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ABSTRACT 


A  full  three  dimensional  computational  study  was  carried  out  using  a  finite- 
volume  based  solver  for  analyzing  the  performance  of  pin-fin  based  micro/meso 
scale  heat  exchangers  with  air  as  the  working  fluid.  A  staggered  arrangement  of 
cylindrical  pin  fins  in  rectangular  channel  geometry  was  used.  Various 
configurations  were  considered  consistent  with  a  parallel  experimental  study 
being  conducted  based  on  a  micro-wind  tunnel  setup.  The  pin/channel  height 
used  was  0.4  mm,  and  the  pin  diameters  varied  from  0.17-0.50  mm  to  give 
hydraulic  diameters  in  the  range  of  0.13-0.78  mm.  This  gave  volumetric  area 
densities  for  the  heat  exchangers  in  the  range  of  5-15  mm2/mm3.  Various  heat 
exchanger  configurations  were  simulated  to  determine  performance 
characteristics  such  as  the  Nusselt  number,  friction  factor,  specific  fluid  friction 
power  and  Mach  number  in  the  Reynolds  number  regime  for  laminar  flows.  In 
addition  a  detailed  numerical  diagnosis  was  carried  out  to  determine  local 
behavior  on  the  pin  surfaces,  end  walls,  etc  to  identify  specific  characteristics 
such  as  regions  of  high  and  low  heat  transfer,  locations  for  possible  shock 
formation,  etc.  The  range  of  results  obtained  would  be  useful  for  future  design  of 
micro  heat  exchangers  for  use  in  small  footprint,  high  heat  flux  dissipation 
applications  like  turbine  blade  and  microelectronic  systems. 
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I.  INTRODUCTION 


A.  BACKGROUND-MOTIVATION 

Today,  great  achievements  and  innovations  are  being  made  in  the 
technological  fields  of  microelectronics  and  micro-electromechanical  (MEMS) 
systems.  The  demands  for  greater  speed,  more  power,  and  less  volume  and 
mass  have  become  more  and  more  urgent  in  most  of  the  forms  and  products  of 
the  science  and  technology.  One  of  the  undesirable  consequences  of  this 
urgency  is  the  operation  in  elevated  temperatures.  Since  the  systems  tend  to 
operate  at  higher  energy  levels,  requirements  are  also  emerging  for  the 
development  of  new  devices  that  can  remove  the  greater  amounts  of  thermal 
energy  and  can  dissipate  the  higher  heat  fluxes.  The  need  for  greater  efficiencies 
and  improved  life  cycles,  which  are  combined  with  less  thermal  stresses, 
accelerated  creep,  and  fatigue  behaviors,  is  growing  too. 

Two  of  the  most  concerned  industries  are  the  microprocessor- 
microelectronics  and  the  gas  turbine  industries.  The  former  has  concentrated  its 
efforts  on  dramatically  reducing  the  size  and  increasing  the  speed  of  its 
attainments.  This  resulted  in  higher  functional  temperatures,  which  create  a 
severe  operational  condition  with  the  significant  effect  of  limiting  the  devices’  life. 
So  the  heat  removal  process  is  critical  and  makes  the  interest  in  micro  heat 
exchangers  essential. 

The  gas  turbine  industry  faces  more  severe  conditions  in  which  micro¬ 
exchanger  technology  could  be  more  applicable,  is  mainly  concerned  with 
increased  performance  which  translates  to  a  higher  inlet  temperature.  Because 
the  turbine’s  reliability  depends  mostly  on  the  mechanical  behavior  of  the  blades, 
their  cooling  due  to  increased  performance  is  extremely  vital.  To  increase  the 
power,  and  thus  the  inlet  temperature,  the  cooling  techniques  must  be  made 
more  effective. 

The  current  methods  that  had  been  used  can  be  divided  into  external  and 
internal.  The  external  method  is  based  on  air  film  cooling  by  which  air  is  injected 
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through  holes  in  the  blade  surface.  The  internal  method  is  based  on  forced 
convection  in  which  air  is  provided  initially  from  the  compressor  bleed  to  the  base 
(mounting)  of  the  blade  and  then  through  an  internal  channel.  A  portion  of  the  air 
can  used  for  the  film  cooling,  and  some  can  be  driven  to  a  pin  array,  which  would 
be  seated  at  the  trailing  edge  which  is  a  suitable  location  (Ref.  1).  The  figure 
below  presents  schematically  the  two  methods. 


Figure  1 . 


Cooled  turbine  airfoil  with  pin  fins  [from  Metzger,  1984] 


In  addition  to  those  two  methods,  other  proposals  include  those  that 
introduce  air  at  the  leading  part  of  the  blade,  ribbed  channels  and  impingement 
cooling  techniques  in  order  to  increase  the  effectiveness  of  the  heat  removal 
process  (Ref.  2).  These  two  techniques  dominated  the  research  for  almost  two 
decades  according  to  a  literature  review.  The  following  figure  shows  two  section 
cuts  of  a  blade  with  the  added  characteristics,  while  the  initial  ones  are  kept 
parallel. 


2 


Tfcrb 

ProMo4«r»  f~~ 


luetba  fulfil 


Figure  2.  Turbine  blade  internal  cooling  passage  [from  Zhang  et  al.,1993] 


A  new  innovative  technique,  which  entails  implementing  an  array  of 
staggered  pins  in  a  blade  shroud  and  then  covering  it  with  a  metal  layer,  was 
recently  introduced.  The  air  after  its  passage  through  the  micro-exchanger 
passages  can  be  driven  outwards  for  external  film  cooling.  Two  manufacturing 
attempts  have  been  reviewed;  one  by  electro-chemical  deposition  (Ref.  3),  the 
other  one  by  LIGA  process  (Ref.  4). 


Figure  3.  SEM  of  the  heat  exchanger  [from  Marques,  2002] 
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Shroud 


Figure  4.  Cross  section  of  heat  exchanger  in  stainless  steel  tube  [from 

Marques,  2002] 


The  above  figures  present  the  configuration  of  the  exchanger  mounted  on 
a  model  blade  and  the  operation  during  which  air  flows  through  the  array  and 
cools  the  hot  blade. 

A  definition  sketch  of  the  chosen  micro-heat  exchanger  array  geometry 
consisting  of  cylindrical  pins  is  shown  in  figure  5.  In  the  geometry,  X  stands  for 
the  streamwise  pitch  between  two  pins,  S  is  the  span-wise  pitch,  H  is  the  pin 
length/height  and  D  is  the  pin  diameter.  This  nomenclature  is  consistent  with  the 
literature. 

Flow 

- > 


Flow 

- > 

Figure  5.  Heat  exchanger  scheme  with  the  common  symbolism  [from  Choo, 

2003] 
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B.  LITERATURE  REVIEW 

Most  of  the  work  in  this  area  so  far  has  been  empirical  and  mainly  in  the 
macro-scale.  Several  researchers  tried  to  understand  the  performance  of 
different  array  configurations,  in  several  conditions,  by  an  examination  of  factors 
like  the  Nusselt  number,  the  heat  transfer  coefficient  and  the  friction  factor.  The 
latter  is  as  significant  as  the  other  heat  transfer  characteristics,  since  in  order  to 
evaluate  a  configuration,  the  pressure  drop  has  to  be  considered.  The  usual  goal 
in  these  studies  is  to  optimize  the  array  to  provide  the  greatest  heat  exchange 
rate  with  the  least  expended  work  for  any  flow  condition. 

One  of  the  early  contributions  came  from  Van  Fossen  (Ref.  5),  who 
studied  the  effect  of  the  presence  of  pins  with  different  H/D  values,  for  a  range  of 
Reynolds  number  from  3,000  -  60,000.  He  concluded  that  the  presence  of  the 
array  returned  higher  heat  transfer  coefficient  values  compared  to  an  empty  duct 
and  they  got  higher  as  the  H/D  was  increased.  His  results  correlated  better  flows 
over  a  Reynolds  number  of  6,000. 

Metzger  et  al.  (Ref.  1,  6)  developed  accurate  correlations  for 
configurations  in  which  the  stream-wise  distance  ratio  was  varied  while  the 
others  were  kept  constant.  Those  correlations  for  the  Nusselt  number  were 
validated  by  Hamilton  (Ref.  7)  with  a  numerical  model.  Also  Metzger  was  one  of 
the  few  who  tried  to  examine  the  effects  of  using  other  types  of  pins,  like  oblong 
shaped,  which  was  also  the  research  subject  of  Arora  (Ref.  8).  Both  Metzger  and 
Arora  concluded  that  the  friction  factor  was  significantly  lower,  but  the  heat 
transfer  rates  did  not  approach  the  cylindrical  pin  values. 

Chen  et  al.  (Ref.  9)  later  examined  more  aerodynamically  shaped  pins  like 
tear  drop  shaped  fins.  They  concluded  that  these  pins  offer  higher  performances, 
and  better  heat  transfer  rates,  with  almost  half  the  friction  factor  of  the  round 
pins.  Those  findings  were  corroborated  by  Hamilton,  who  found  that  the  optimum 
airfoil  shaped  pin  array  outperforms  a  similar  cylindrical  pin  array,  which  needs 
triple  the  specific  energy  loss  to  produce  analogous  heat  transfer  rates. 
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There  has  been  considerable  deliberation  among  the  researchers  with 
regard  to  the  heat  fluxes  from  the  pins  compared  to  the  endwalls.  Chyu  et  al. 
(Ref.  10)  show  comparable  fluxes  between  these  two  parts  of  the  exchanger; 
unlike  Metzger  and  Van  Fossen,  who  showed  that  the  pins  generated  50  percent 
higher  heat  transfer  rates  than  the  end-walls.  On  the  other  hand,  Al  Dabagh  et  al. 
(Ref.  1 1 ),  found  that  end-walls  show  heat  transfer  rates  35  percent  higher  than 
the  staggered  pins.  Finally,  Chyu  et  al.  (Ref.  12)  in  alter  study  found  that  the  pin 
surfaces  offer  around  20  percent  more  in  heat  transfer  than  the  endwalls. 
Furthermore  their  results  for  the  Nusselt  number  are  very  close  to  those  reported 
by  Metzger  (Ref.  6) 

Also  very  interesting  is  the  approach  of  Tahat  et  al.  (Ref.  13,  14),  who  tried 
to  correlate  the  Nusselt  number  with  the  Reynolds  number,  including  the 
geometry  configuration  of  the  tested  arrays  by  using  the  ratios  of  the  stream  and 
the  span-wise  distance  of  the  pins  over  the  total  length  and  width  of  the  array. 
They  used  a  specific  test  bed  with  a  varying  ceiling  height  to  test  different  pin 
bank  arrangements  to  find  the  optimum  spacing  for  the  pins. 

In  1997,  Li  et  al.  (Ref.  15)  conducted  an  experimental  work  with  elliptical 
pin  arrays  and  compared  their  data  with  Metzger’s  findings.  They  used  the  same 
dimension  pins  in  two  different  arrays,  but  with  equal  stream  and  span-wise 
distance  between  them,  for  a  range  of  the  Reynolds  number  800  -  9,000.  They 
found  better  heat  transfer  rates  with  considerable  lower  pressure  losses.  The 
more  packed  his  configuration  was,  the  greater  the  performance  with  resulting 
higher  Nusselt  numbers. 

One  of  the  most  recent,  and  probably  the  first  in  micro  dimensions 
experimental  work  is  by  Marques  et  al  (Ref.  16),  who  tested  the  configuration  of 
X/D=S/D=2.5  and  H/D=1  for  Reynolds  number  from  4,500  to  19,200.  The  same 
configuration  was  also  tested  by  Metzger  and  Chyu  in  the  past,  but  in  the  macro¬ 
scale.  Marques  et  al.  developed  a  correlation  for  the  Nusselt  number  for  the 
range  of  Reynolds  noted  above,  which  gives  more  optimistic  results  than  the 
other  two  researchers.  A  possible  explanation  is  that  their  theoretical  model  uses 
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a  smaller  heat  exchange  area  in  the  array  than  the  exchanger  footprint  area. 
They  then  developed  a  theoretical  model  that  can  compute  the  performance  of 
such  a  heat  exchanger  when  it  is  used  under  a  shroud  in  a  turbine  blade  to 
cooling  it. 

Finally,  at  the  Naval  Postgraduate  School,  Monterey,  California,  in  2002- 
04  an  experimental  apparatus  was  built  for  testing  a  set  of  compact  heat 
exchangers  by  Roussakies  (Ref.  17)  and  Bosserman  (Ref.  18).  The  arrays  are 
the  same  as  those  that  will  be  studied  in  this  work.  Because  some  design  issues 
of  the  testing  bed  facility  arose,  the  results  are  not  as  reliable  as  had  been 
expected  which  a  good  reason  for  also  checking  them  against  a  numerical 
model. 

On  the  other  hand,  little  numerical  modeling  work  had  been  done  until 
recently  in  comparison  with  empirical.  Recent  advances  in  CFD  tools  have  made 
this  a  very  promising  avenue  of  research.  Two  important  factors  in  the  push 
toward  that  direction  are  the  great  manufacturing  and  testing  difficulties  and  the 
high  cost  that  these  studies  demand. 

Donahoo  et  al.  (Ref.  19)  developed  a  two  dimensional  model  of  staggered 
pin  arrays,  by  using  a  general  purpose  viscous  solver,  in  order  to  optimize  their 
configurations  in  axial  pitch  (X/D)  spacing  for  turbulent  Reynolds  numbers.  They 
found  that  the  maximum  heat  transfer  occurred  between  the  fourth  and  sixth  row 
for  the  tested  arrays.  Finally  they  recognized  that  a  3  D  model  is  a  future 
necessity  since  it  was  not  possible  to  capture  the  pin-endwall  interactions  with 
the  2  D  model  and  compare  their  participation  to  the  whole  heat  transfer  of  the 
array  or  to  study  the  effect  of  the  H/D  ratio. 

Flamilton  (Ref.  7)  used  the  finite  element  commercial  code  “ANSYS”  to 
develop  a  numerical  3  D  model  to  study  several  configurations  of  staggered  pin- 
fin  heat  exchangers,  but  still  in  the  macro-scale.  Initially,  he  validated  his  model 
with  prior  experimental  results  from  Arora  et  al.  and  Metzger  et  al.  He  then 
proceeded  to  examine  of  the  performance  of  arrays  with  cylindrical  pins.  He 
developed  correlations  for  the  Nusselt  number  for  all  his  configurations,  and  also 
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found  that  pin  surfaces  had  from  50-200  percent  higher  heat  transfer  coefficients 
than  the  walls.  Achieving  small  increases  of  it  requires  disproportional  increases 
in  friction  power  losses.  He  concluded  that  variations  in  axial  pitch  did  not  have 
as  great  an  effect  in  performance  as  those  the  in  H/D  ratio.  He  also  studied 
airfoil-shaped  pin  arrays,  from  which  he  concluded  that  they  are  better 
performers  than  corresponding  dimensional  cylindrical  pin  arrays  since,  to  attain 
comparable  heat  flux  removal,  they  require  only  a  third  of  the  energy  as 
compared  to  the  latter. 

Choo  (Ref.  20)  carried  out  a  laminar  flow  numerical  analysis  in  2003.  He 
simulated  exchangers  performing  in  micro-scale  in  the  laminar  regime  with 
Reynolds  numbers  in  the  range  of  100-1,000  to  study  the  effect  of  the 
geometrical  ratios  of  the  exchangers  (X/D,  S/D,  H/D).  He  stated  that  for  his 
geometries  and  simulations  the  maximum  heat  transfer  occurred  around  the 
second  row  of  the  models.  Also,  after  taking  into  account  all  three  effects,  he 
suggested  that  the  optimum  configuration  of  his  matrix  was  X/D=1.25  and 
S/D=H/D=3. 

To  minimize  the  amount  of  time  and  cost  that  are  needed  in  an  empirical 
project,  use  of  a  computational  study  for  pre-evaluation  and  optimization  of 
several  designs  is  absolutely  essential.  Another  characteristic  worth  mentioning 
is  the  tremendous  flexibility  that  a  computational  study  can  offer  as  compared  to 
empirical  methods.  Finally,  the  attributes  of  the  theory  that  were  developed  for 
the  macro-scale  studies  can  also  be  directly  applied  to  micro-scale,  since  they 
are  based  on  the  same  characteristic  qualitative  factors. 

C.  OBJECTIVES 

The  main  objective  of  this  study  is  to  determine  the  performance  and 
characteristics  of  flow  and  heat  transfer  behavior  of  a  set  of  micro  heat 
exchangers  of  specific  configurations.  A  3-D  model  will  be  developed  that  will  be 
validated  by  theoretical  and  historical  experimental  results  in  order  to  use  it  to 
examine  the  behavior  of  the  different  staggered  cylindrical  pin  arrays  in  a  laminar 
compressible  flow  regime. 
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The  parameters  that  will  be  inspected  are  the  Nusselt  number,  the  heat 
transfer  coefficient,  the  Mach  number,  the  array  pressure  drop,  the  friction  factor, 
and  the  specific  fluid  friction  power  for  different  Reynolds  numbers  .The 
theoretical  definitions  for  the  above  parameters  will  be  provided  in  the  following 
chapter.  The  results  will  be  used  to  define  suitable  regions-conditions  of 
operation  for  every  heat  exchanger’s  geometry  and  will  be  compared  to  findings 
of  experimental  tests  already  completed. 


D.  METHOLOGY -ORGANIZATION 

Several  three-dimensional  models  were  constructed  according  to  our  test 
matrix,  using  the  commercial,  computational  fluid  dynamics  package  CFD-ACE+. 
The  construction  was  carried  out  with  the  preprocessor  of  the  package  CFD- 
GEOM.  Then  the  models  were  transferred  to  the  processor  CFD-ACEU  for 
executing  the  simulation.  Finally,  the  postprocessor  of  the  package  CFD-VIEW 
was  used  to  study,  analyze,  and  manipulate  the  results.  The  details  of  the 
modeling  will  be  discussed  in  the  next  chapter.  In  the  third  chapter  we  provided 
the  validation  and  corroboration  with  the  theoretical  and  historical  results  and  a 
more  extensive  heat  exchange  examination  of  the  observations.  In  the  fourth 
chapter  the  results  of  the  simulations  of  the  tested  configurations  are  discussed. 
The  final  chapter  draws  certain  conclusions  and  provides  some 
recommendations. 
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II.  FINITE  VOLUME  MODELING 


A.  MODELING  GENERAL  DESCRIPTION 

The  CFD-ACE+  package  solver  was  used  in  this  study.  It  is  a  solver  that 
uses  a  finite  volume  pressure  correction  based  method  to  solve  the  conservation 
equations  of  mass,  momentum  and  energy.  The  software  provides  the  user  with 
three  graphical  user  interfaces  (GUIs):  GEOM  for  the  geometries-models 
creation,  the  ACEU  that  is  the  actual  solver  of  the  discretized  equations,  and  the 
postprocessor  VIEW  for  viewing  and  post-processing  the  results.  CFD-ACE+  has 
been  very  popular  for  this  kind  of  problem,  such  as  fluid  flow  fields  in  conjugation 
with  heat  transfer,  because  it  provides  reliable  results  for  a  variety  of  problems. 
This  independence  makes  this  solver  more  powerful  since  it  provides  very  similar 
results  for  several  grid  refinements  within  good  accuracy. 


B.  TEST  MATRIX 

The  test  matrix  of  this  study  consists  of  ten  different  configurations.  The 
first  eight  configurations  and  the  tenth  have  been  tested  here  at  the  Naval 
Postgraduate  School  by  testing  Bosserman  and  Roussakies;  but  because  of 
some  construction  constraints  and  limitations  of  the  utilized  bed,  their  data  is  still 
preliminary  and  is  being  verified.  All  of  the  configurations  have  the  same  height 
due  to  the  specific  height  of  the  flow  channel  inside  the  designed  apparatus 
(H=405pm).  The  first  four  have  a  pin  diameter  of  500pm  (D=500pm),  while  the 
next  set  of  four  has  a  pin  diameter  of  166.7pm  (D=1 66.7pm).  The  tenth  one  has 
a  pin  diameter  of  151.5pm;  it  was  Summers  (Ref.  21)  optimum  laminar 
configuration  in  the  macro-scale  level.  The  two  sets  of  four  vary  under  the  same 
ratios  for  the  stream-wise  and  the  span-wise  pin  spacing.  This  can  be  observed 
in  the  following  table,  which  presents  the  tested  matrix,  with  some  geometric 
characteristics  that  will  be  explained  thoroughly  later  in  this  chapter. 
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Table  1 .  Test  matrix  for  the  heat  exchanger  configurations 


# 

X/D 

S/D 

H/D 

D  (pm) 

aH 

Dh(pm) 

V0pen(m  ) 

AWetted(m  ) 

HX#1 

3 

3 

0.81 

500 

2.11 

701 

8.32E-10 

4.74E-06 

HX#2 

1.25 

3 

0.81 

500 

2.26 

567 

3.00E-10 

2.12E-06 

HX#3 

1.25 

1.25 

0.81 

500 

2.62 

307 

7.87E-1 1 

1.02E-06 

HX#4 

3 

1.25 

0.81 

500 

2.26 

567 

3.00E-10 

2.12E-06 

HX#5 

3 

3 

2.43 

167 

2.67 

553 

9.24E-1 1 

6.68E-07 

HX#6 

1.25 

3 

2.43 

167 

3.62 

354 

3.34E-1 1 

3.77E-07 

HX#7 

1.25 

1.25 

2.43 

167 

5.88 

137 

8.74E-12 

2.55E-07 

HX#8 

3 

1.25 

2.43 

167 

3.62 

354 

3.34E-1 1 

3.77E-07 

HX#9 

2.5 

2.5 

1 

500 

2.25 

777 

6.83E-10 

3.52E-06 

HX#10 

5 

5 

2.67 

152 

2.27 

690 

2.25E-10 

1.30E-06 

The 

ninth 

configuration 

is  the 

one 

used  most  by  the 

researchers, 

especially 

in  the 

macro-scale 

and  it  is 

i  the 

one  that 

is  going  to 

be  used  for 

validations  and  corroborations,  along  with  the  seventh,  which  has  the  largest 
magnification  factor  (aH)  and  is  the  “densest”  of  all  of  them.  The  area  density  of 
an  array  unit  cell  is  defined  as  the  ratio  of  the  wetted  (heat  transfer)  surface  area 
over  the  total  volume.  The  magnification  factor  that  is  provided  in  the  table  is 
defined  as  the  product  of  the  area  density  and  the  pin  height.  For  an  exchanger 
without  pins  that  is  a  rectangular  cross  sectional  duct,  the  aH  becomes  equal  to 
2.  All  the  above  table  values  have  been  computed  based  on  a  unit  cell. 


Table  2.  Unit  cell  dimensions  for  the  tested  heat  exchangers 


# 

X/D 

S/D 

H/D 

D  (pm) 

S  (pm) 

X  (pm) 

H(pm) 

HX#1 

3 

3 

0.81 

500 

1500 

1500 

405 

HX#2 

1.25 

3 

0.81 

500 

1500 

625 

405 

HX#3 

1.25 

1.25 

0.81 

500 

625 

625 

405 

HX#4 

3 

1.25 

0.81 

500 

625 

1500 

405 

HX#5 

3 

3 

2.43 

167 

500 

500 

405 

12 


HX#6 

1.25 

3 

2.43 

167 

500 

208 

405 

HX#7 

1.25 

1.25 

2.43 

167 

208 

208 

405 

HX#8 

3 

1.25 

2.43 

167 

208 

500 

405 

HX#9 

2.5 

2.5 

1 

500 

1250 

1250 

500 

HX#10 

5 

5 

2.67 

152 

758 

758 

405 

The  dimensions  of  the  unit  cells  for  all  the  configurations  are  shown  in  the 
Table  2.  The  heat  exchangers  have  been  divided  into  two  groups  and  they  will  be 
examined  according  to  their  pin  diameter. 


C.  BOUNDARY  CONDITIONS 

Several  assumptions  had  to  be  made  for  the  purpose  of  reducing  the 
computational  requirements.  Because  the  geometry  and  scale  of  the  micro  heat 
exchanger  it  usually  contains  a  large  number  of  pins,  and  the  contribution  of  pins 
as  compared  to  the  endwalls  in  heat  transfer  will  be  sufficiently  greater.  Micro 
heat  exchangers  have  a  very  large  number  of  pins  in  the  span-wise  direction, 
and  the  contribution  of  the  side  walls  in  the  heat  exchange  process  could  be 
relatively  neglected.  Instead  Hamilton  found  that  it  is  sufficient  to  consider  only  2 
spanwise  rows  and  introduce  symmetry  sidewall  boundary  conditions  to  account 
for  the  ignored  rows.  A  symmetric  boundary  was  also  assumed  at  the  Hx 
midplane  in  order  to  simulate  only  half  of  the  geometry.  The  symmetric  planes 
have  the  zero  cross-over  velocity  components,  and  from  the  thermal  perspective 
is  are  adiabatic  in  nature.  Furthermore,  the  wall  boundary  condition  was  used  for 
the  pin  surfaces  and  the  lower  endwall  where  the  no  slip  condition  was  applied. 

From  the  heat  transfer  point  of  view,  the  pins  and  the  lower  endwall  were 
assumed  to  be  isothermal.  A  justification  calculation  for  the  isothermal  character 
of  the  pins  can  be  found  in  Hamilton  (Ref.  7).  An  illustration  of  the  model  is 
presented  in  figure  6. 
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Heated  pin 
Surface 


Isothermal 

Heating 


Figure  6.  Schematic  of  the  numerical  model 


To  validate  the  models  with  the  theoretical  and  historical  experimental 
results,  an  entry  and  an  exit  straight  duct  section  was  used.  The  entry  duct  is 
used  to  provide  a  fully  developed  velocity  profile  at  the  flow  entrance  to  the 
exchanger;  while  the  exit  duct  will  provide  well  mixed  conditions  for  the  flow  so  as 
to  make  the  exit  temperature  calculations  more  effective  and  robust. 
Subsequently,  we  used  the  concept  of  the  entry  and  exit  ducts  in  all  the  array 
configurations  in  order  to  be  consistent.  In  both  ducts  we  applied  the  same 
boundary  conditions  for  the  surrounding  planes.  The  sidewalls  and  the  upper 
plane  were  assumed  to  be  adiabatic  symmetry  planes,  while  the  bottom  was 
assumed  to  be  an  adiabatic  wall.  Furthermore,  we  assumed  an  inlet  with  a  fixed 
velocity  that  fixes  the  mass  flow  rate  and  thus  the  Reynolds  number  of  the 
simulation  and  a  fixed  pressure  outlet  boundary.  A  representation  of  the  above 
conditions  is  provided  in  the  figure  7. 
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Heated  pin  surface 
at  T=320K 


Heated  endwall 
at  T=320K 


Air  inlet 
T=300K 


Adiabatic  exit  endwall 


Entry  duct 


Adiabatic  Symmetry  Condition,  Uy=Uz=0 


Adiabatic  entry  endwall 


Figure  7.  Boundary  conditions  representation 


D.  GRID  GENERATION 

Since  we  decided  to  use  the  “GUI”  of  the  GEOM  and  not  the  PYTHON 
scripting  feature,  we  had  to  create  several  geometrical  models  independently. 
First,  the  2-D  model  of  a  half  unit  cell  had  to  be  constructed,  since  we  intended  to 
use  an  unstructured  triangular  meshing.  Thus  the  creation  of  a  closed  loop  was 
unavoidable.  The  next  step  was  to  duplicate  that  half  unit  cell  at  the  proper 
position,  according  to  every  array  configuration.  We  did  a  very  careful  meshing 
that  was  finer  near  the  pin’s  perimeter  and  coarser  towards  the  internal  part,  to 
capture  the  large  heat  and  flow  gradients  near  the  walls.  In  most  cases  we  used 
a  30  degrees  curvature  criterion  and  a  growing  factor  from  1 .05  to  1.10.  Then  we 
linearly  extruded  the  whole  geometry  toward  the  third  direction  by  using  a  power 
law  with  the  purpose  again  to  create  layers  of  triangular  prisms  that  will  be  finer 
near  the  heated  endwall  where  the  gradients  are  also  expected  to  be  higher.  We 
carefully  examined  all  the  3-D  centroid  /  face  angles  to  ensure  that  were  at  least 
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25  degrees  and  also  that  no  negative  volume  cells  had  been  created.  Figure  8 
provides  a  fingerprint  illustration  of  the  grid  around  a  unit  cell  of  an  array. 


Figure  8.  Sample  unit  cell’s  meshing  in  2-D 


Another  issue  was  the  entry  length  in  the  ducts.  Since  we  had  to  construct 
several  configurations  and  run  them  in  different  Reynolds  numbers  we  had  to 
figure  out  a  way  to  reduce  the  model  generation  work  and  eliminate  the  need  for 
different  ducts  for  every  simulation.  This  was  achieved  by  using  a  feature  of  the 
solver  ACE  where  we  can  simply  input  the  velocity  profile  for  the  inlet  boundary 
to  ensure  fully  developed  conditions.  Therefore,  we  used  a  constant  length  duct 
for  all  the  models  that  were  generated,  with  lengths  of  three  and  five  hydraulic 
diameters  for  the  entry  and  exit  sections  respectively.  Since  the  height  of  all  the 
models  is  the  same,  the  duct  hydraulic  diameter  for  all  the  models  is  constant 
and  equal  to  twice  the  height  for  an  infinite  span-wise  model. 
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Because  we  used  structured  grids  for  meshing  the  two  ducts,  the  whole 
model  had  a  hybrid  character,  and  we  tried  not  to  exceed  a  block  aspect  ratio  of 
five  in  the  cell  extrusion  process.  This  would  serve  the  connection  regions 
between  the  two  grids  for  having  similar  cell  dimensions  to  ensure  a  smooth 
transition  from  one  volume  to  the  next. 


Figure  9.  Sample  model’s  grid 


Figure  8  illustrates  a  sample  mesh  of  a  model.  After  the  generation  of  any 
geometry  we  had  to  scale  it  in  micrometers  in  order  to  import  that  into  the  solver 
GUI. 

E.  SOLUTION  METHOD  -  SOLVER  CONTROLS 

As  mentioned  previously  the  simulations  were  executed  in  the  laminar 
regime.  Hamilton  (Ref.  7)  in  his  dissertation  showed  that  the  transition  typically 
occurs  at  a  Reynolds  number  1000-1500.  Within  that  region,  he  achieved  an 
overlap  in  a  figure  of  the  laminar  and  turbulent  models  that  he  used.  Zhukauskas 
(Ref.  23)  makes  a  similar  statement  in  his  research  review.  The  Reynolds 
number  on  which  our  runs  are  based  extended  from  100  to  1 ,000. 
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In  the  CFD-ACE+  solver  the  flow  and  heat  transfer  modules  were 
selected  in  order  to  obtain  the  flow  (velocity  and  pressure)  and  temperature 
distributions  for  the  fluid  volume  elements.  To  formulate  an  equation  for  the 
unknown  variable  pressure,  the  continuity  equation  is  solved  by  adopting  the 
SIMPLEC  algorithm  in  CFD  -  ACE  to  start  the  iteration  process.  For  more 
information,  the  reader  is  referred  to  the  user  manual  (Ref.  24)  or  the  book  by 
Versteeg  and  Malalasekera  (Ref.  25). 

In  regard  to  the  volume  conditions,  because  we  are  going  to  solve  the 
compressible  form  of  the  numerical  models,  the  density  is  defined  from  the  ideal 
gas  law,  which  is  the  sixth  equation  in  the  governing  equation  set.  Furthermore, 
the  air  viscosity  was  obtained  by  using  Sutherland’s  law  for.  Also,  the  air 
pressure  coefficient  Cp  was  held  constant,  since  its  variation  in  the  region  of  300 
-  320  K  is  negligible.  For  the  thermal  conductivity,  which  varies  more  than  the 
pressure  coefficient,  we  used  the  constant  Prandtl  number  law  which  is  a 
variation  of  k  according  to  the  definition  of  the  Prandtl  number  which  is  held 
constant  at  a  value  of  0.707. 


Pr  = 


k  = 


(1) 


The  next  step  in  the  solver  tabs  is  to  input  the  boundary  conditions.  The 
velocities  for  the  wall  boundaries,  pins  and  endwall  were  kept  at  zero  and  the 
temperature  was  set  to  320  K.  The  inlet  temperature  was  set  to  300  K,  and  for 
the  flow  part  we  input  a  fully  developed  X-velocity  profile  in  the  Z  direction.  The 
profile  for  every  Reynolds  number  and  for  every  array  was  calculated  by  using  an 
ECXEL  spreadsheet;  which  was  then  utilized  by  the  solver  to  fix  the  mass  flow 
rate  which  subsequently  derives  the  Reynolds  number  for  this  simulation. 

For  the  flow  initial  conditions  we  input  a  velocity  that  was  based  on  the 
duct  cross  sectional  area  and  was  derived  by  dividing  the  specific  mass  flow  rate 
with  by  duct  area  and  film  density,  while  the  initial  temperature  was  set  to  300  K. 

Another  crucial  control  is  the  scheme  that  we  had  to  choose  in  order  to 


discretize  the  advection  term,  especially  when  this  term  dominates  other  terms  in 
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the  governing  equations.  The  monotone  streamline  upwind  (MSU)  scheme  is 
generally  quite  robust  and  produces  stable  diagonally  dominant  matrixes,  but  is 
first  order  accurate.  It  is  also  the  scheme  recommended  in  the  ACE+  manual, 
especially  for  complex  geometries. 

Trial  runs  were  carried  out  using  other  second  order  schemes  like  the 
central,  or  the  upwind  limiter,  always  blended  with  the  MSU  to  prevent 
divergence.  The  same  scheme  is  also  suggested  by  the  solver  manual 
documentation.  We  tried  to  overcome  the  choice  of  the  first  order  scheme  by 
increasing  the  convergence  criteria  for  the  residuals  average  and  by  using  more 
refined  models  for  the  simulation.  As  an  example  of  the  deviation  between  the 
schemes  for  HX#7,  and  for  the  case  of  Re=100,  the  percentage  difference 
between  the  upwind  and  the  central,  as  far  as  it  concerns  the  heat  addition  to  the 
fluid  was  0.2  percent  and  between  the  upwind  limiter  was  0.3  percent. 

As  far  as  it  concerns  the  relaxation  criteria,  in  most  cases  the  solver 
default  values  worked  well,  except  for  the  higher  Reynolds  number  simulations, 
near  1,000,  where  we  had  to  modify  them  a  little  by  increasing  the  inertial 
relaxation  values  up  to  0.3  -  0.4  and  by  decreasing  the  linear  ones  down  to  0.9  in 
order  to  provide  more  stability  to  the  solution.  This  slowed  down  the  convergence 
process  but  made  it  achievable. 

F.  NUMBER  OF  PIN  ROWS 

Another  issue  that  we  had  to  resolve  was  the  number  of  pin  rows  that  we 
were  going  to  use  for  every  array  that  we  were  going  to  model.  Because  of  the 
great  variation  of  the  configurations,  a  great  variation  in  their  performances  was 
expected.  For  every  configuration,  and  especially  for  the  lower  Reynolds  number 
simulations,  the  bulk  outlet  temperature  was  expected  to  be  higher.  It  began  to 
reduce  as  the  Reynolds  number  was  increased,  which  agrees  with  the  theory,  in 
Kays  and  Crawford  “Convective  Heat  and  Mass  Transfer”  (Ref.  26).  So  for  every 
array  we  had  to  use  only  so  many  rows  of  pins  that  would  yield  an  outlet  bulk 
temperature  no  more  than  319  K  to  319.5  K,  or  0.5  to  1  K  lower  than  the 
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prescribed  wall  and  pins  temperature  of  320  K,  in  order  to  avoid  errors  due  to 
loss  of  resolution  in  the  numerical  calculations. 

Therefore,  for  every  configuration  we  generated  a  geometry  without  ducts 
and  with  ten  rows  of  pins,  and  we  ran  that  for  Reynolds  of  100.  Then  by  utilizing 
CFD  -  VIEW,  we  constructed  a  carpet  plot  of  the  temperature  for  the  upper 
symmetry  plane,  illustrated  in  figure  10,  to  get  an  initial  rough  estimate  in  which 
row  we  achieved  a  temperature  around  319  K. 


Figure  10.  Temperature  profile  for  two  pin  arrays  with  ten  rows 


Then  for  higher  accuracy  we  went  to  this  row  of  fluid  blocks,  made  a 
stream  wise  direction  (X)  cut  and  integrated  the  results  on  that  plane.  This  was 
done,  to  find  an  accurate  bulk  temperature  at  that  exact  X  coordinate.  By  moving 
the  cut  in  the  selected  row  blocks  we  could  check  the  temperature  continuously 
in  every  position. 

Then  we  went  back  to  GEOM  to  generate  a  final  version  of  the  model  with 
the  desired  number  of  pin  rows  and  with  the  entry  and  exit  duct  as  it  were 
described  above  for  the  specific  runs. 
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Table  3.  Number  of  pin  rows  used  per  configuration 


# 

X/D 

S/D 

H/D 

Nx 

Ns 

L(m) 

W(m) 

HX#1 

3 

3 

0.81 

6 

2 

0.0090 

0.0030 

HX#2 

1.25 

3 

0.81 

10 

2 

0.0063 

0.0030 

HX#3 

1.25 

1.25 

0.81 

4 

2 

0.0025 

0.0013 

HX#4 

3 

1.25 

0.81 

4 

2 

0.0060 

0.0013 

HX#5 

3 

3 

2.43 

8 

2 

0.0040 

0.0010 

HX#6 

1.25 

3 

2.43 

10 

2 

0.0021 

0.0010 

HX#7 

1.25 

1.25 

2.43 

4 

2 

0.0008 

0.0004 

HX#8 

3 

1.25 

2.43 

4 

2 

0.0020 

0.0004 

HX#9 

2.5 

2.5 

1 

6 

2 

0.0075 

0.0025 

HX#10 

5 

5 

2.67 

8 

2 

0.0061 

0.0015 

Table  3  provides  the  specific  number  of  rows  that  every  heat  exchanger 
configuration  used.  All  the  simulations  used  two  unit  cells  in  the  span-wise 
direction.  In  the  above  table  the  dimensions  of  the  configurations  for  all  the  heat 
exchangers  are  also  presented. 


G.  THEORETICAL  APPROACH 

1.  Characteristic  Length 

The  most  appropriate  characteristic  length  that  is  the  hydraulic  diameter  of 
this  Hx  array  is  obtained  from  flow  theory  in  porous  media,  and  is  defined  as 
follows  in  terms  of  the  open  volume  that  is  available  for  the  fluid  flow  in  the  array 
and  the  wetted  heat  transfer  surface  area: 

4*  V 

Dh=_A~~^  (2) 

''wetted 

It  is  the  most  truly  dimensionless  quantity  that  represents  consistently  the 
tortuous  flow  paths  inside  the  exchanger.  It  is  the  same  characteristic  length  that 
was  proposed  by  Van  Fossen  (Ref.  5)  and  that  is  used  by  many  more  recent 
experimental  and  numerical  researchers.  The  hydraulic  diameter  remains  the 
same  whether  its  calculation  are  made  based  on  the  unit  cell  or  based  on  the 
whole  array. 
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2.  Entry  -  Exit  Length  Considerations 

According  to  White  (Ref.  22),  for  every  Reynolds  number  the  relation  that 
gives  the  entrance  length  for  a  fully  developed  flow  is  the  following: 

Le  «0.06*Dh  *ReDh  (3) 

Because  it  was  not  possible  to  construct  a  different  entry  duct  for  every 
simulation  and  because  the  existence  of  ducts  usually  smoothes  the  flow 
conditions,  we  decided  to  impose  a  fully  developed  velocity  profile  and  to  use  the 
constant  length  specified  above  for  all  the  flow  conditions. 

For  the  exit  duct  length,  which  is  more  crucial  because  it  is  necessary  for 
ensuring  a  well  mixed  flow  that  is  starting  to  redevelop  again  and  because  it  will 
help  in  a  more  accurate  integration  for  computing  the  predicted  outlet  bulk 
temperature,  we  used  a  little  longer  section  than  the  entry  one,  equal  to  five 
hydraulic  diameters. 

3.  Inlet  Velocity  Considerations 

According  to  White  (Ref.  22),  the  relation  that  gives  the  velocity  profile  for 
a  Poiseuille  flow  between  parallel  plates  is 

u  =  (4) 

where  umax  is  the  maximum  velocity  at  the  centerline  of  the  duct,  which  is 
equal  to 

Un,„=fVi„(5) 

where  Vin  is  the  average  inlet  duct  velocity.  Consequently  by  using  the 
equation  (4)  we  created  a  six  -  point  profile  along  the  z-axis,  which  then  was 
used  at  the  definition  of  the  inlet  boundary  condition  in  the  solver. 

4.  Reynolds  Number 

The  Reynolds  number  that  we  are  going  to  use  is  based  at  the  previously 
defined  hydraulic  diameter  and  in  a  characteristic  velocity,  the  average  array 
velocity 
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(6) 


Re  = 


pUDh 

M 


where 


(7) 


and  A  is  the  average  array  cross  sectional  area,  which  is  equal  to 

-  V 

A  _  _^erL  (8) 

Also,  according  to  the  above  relations,  the  Reynolds  number  can  be 
expressed  associated  to  the  mass  flow  rate.  This  relation  will  be  used  after  every 
simulation  to  calculate  the  exact  Reynolds  number  from  the  “mass  flow 
summary”  from  the  output  file  of  ACE. 


Re  = 


m£h 

pA 


(9) 


H.  COMPARED  PARAMETERS  ANALYSIS 

1.  Nusselt  Number 

The  Nusselt  number  is  a  dimensionless  number  that  provides  the  ratio  of 
the  convective  heat  transfer  over  a  surface  that  would  occur  by  fluid  motion  to  the 
corresponding  conductive  heat  transfer  and  is  a  measure  of  the  heat  transfer 
ability.  In  our  case  more  meaning  makes  the  average  Nusselt  number  which  is 
given  by  the  following  relation: 

Nu= (10) 
k 

where  harray  is  the  average  array  heat  transfer  coefficient  and  can  be 

calculated  by  using  the  thermal  energy  difference  between  the  flow  inlet  and  the 
outlet,  which  is  provided  by  the  heat  transfer  summary  in  the  output  file  of  ACE  at 
the  end  of  the  simulation,  and  by  using  the  log  mean  temperature  difference  as 
below: 
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AQ 


(11) 


h 


array 


A  AT 

''wetted 4-1  1  Im 


For  computing  the  log  mean  temperature  difference  we  will  use  the 
definition  of  the  outlet  bulk  fluid  temperature,  which,  according  to  Kays  and 
Crawford  in  “Convective  Heat  and  Mass  Transfer”  (Ref.  26)  is  provided  from 


ml 


bulk, out 


(12) 


In  the  above  definition  the  Cp  is  not  included  because  we  had  assumed 
already  that  it  is  constant  and  does  not  influence  the  relation.  The  outlet  bulk  fluid 
temperature  was  calculated  in  the  VIEW  program  by  integrating  at  the  outlet  of 
the  adiabatic  exit  duct  with  respect  to  the  duct  cross  sectional  area  by  using  the 
expression  calculator  feature.  It  is  the  temperature  that  characterizes  the  average 
thermal  energy  of  the  flow,  since  it  is  based  at  the  mass  flow  rate. 

The  log  mean  temperature  difference  is  given  by  the  formula  below; 


A-r  _  ("l"wall  "l"bulk,in  )-(Tw 
nl|m  fT  _  T 


where  the  inlet  bulk  fluid  temperature  was  constant  and  equal  to  300K. 

2.  Effective  Heat  Transfer  Coefficient 

The  effective  heat  transfer  coefficient  is  the  corresponding  heat  transfer 
coefficient  in  an  empty  (without  pins)  open  channel  had  to  provide  in  order  to 
transfer  the  same  amount  of  thermal  energy  like  the  one  with  the  pins. 

u  u  aH  . .  .. 

heff  =harray^-  (I4) 

It  is  an  average  quantity  of  the  whole  array  since  it  is  based  at  the  average 
array  heat  transfer  coefficient,  and  it  becomes  equal  to  it  for  the  open  channel 
because  aH=2. 
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3.  Friction  Factor 

The  dimensionless  friction  factor  of  the  pin  array  that  was  used  was 
defined  as  follows: 

f  =  AParrayDh  (15) 

— pU2L 
2 

where  AParray  was  the  total  pressure  difference  between  the  inlet  and 
outlet  sections  of  the  array.  It  was  computed  from  the  VIEW  post  processor  by 
integrating  at  the  sections  with  the  expression  calculator  tool  in  order  to  find  an 
average  value  of  the  inlet  and  outlet  section.  L  is  the  total  stream-wise  length  of 
the  array.  Also,  the  density  that  was  used  was  the  one  based  on  the  standard 
atmospheric  conditions  but  calculated  in  the  film  temperature,  in  order  to  have 
the  same  measure  of  comparison  between  the  different  heat  exchanger 
configurations. 

4.  Specific  Fluid  Friction  Power 

Kays  and  London  (Ref.  27)  define  a  quantity  that  can  express  the 
performance  of  a  heat  exchanger,  the  specific  fluid  friction  power,  which 
characterizes  the  work  that  must  be  expended  to  overcome  the  viscous  effects 
per  unit  wetted  area. 

E  =  mAF>array  (16) 

Pfilm^  wetted 

5.  Effectiveness  -  NTU  Relations 

According  to  Incropera  and  De  Witt  in  “Introduction  to  Heat  Transfer”  (Ref. 
28)  the  effectiveness  of  a  heat  exchanger  is  defined  as  the  ratio  of  the  actual 
heat  transfer  rate  to  the  maximum  possible  heat  transfer  rate  and,  in  our  case,  is 
defined  as: 


£  = 


"^"bulk.out  "^"bulk.in 
~l"wall  "^"bulk.in 


(17) 
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The  same  source  shows  that,  for  all  the  exchangers,  the  effectiveness 
associated  with  the  number  of  transfer  units  (NTU)  is  given  by  the  following 
relation: 


s  =  l-exp(-NTU)  (18) 
The  NTU  generally  is  defined  as 

I  )A 

NTU  =  — —  (19) 


but  in  our  case  as 


1 

UA 


array 


1_ 

A 


wetted 


(20)  and  Cmin  =  rhCp  (21) 


6.  Maximum  Mach  Number 

The  maximum  Mach  number  was  extracted  with  the  post  processor  by 
examining  the  fluid  volume  inside  the  array  and  by  again  using  the  expression 
calculator  and  the  max  built  in  function.  In  order  to  define  suitable  operating 
regions  for  all  the  exchangers,  the  maximum  Mach  number  will  be  plotted  versus 
the  corresponding  Reynolds  number  for  comparison. 
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III.  VALIDATION  -  CORABORATION  AND  TRANSITION 

EXAMINATION 


A.  GRID  INDEPENDENCE  STUDY 

It  was  anticipated  that  the  simulations  of  these  micro-scale  heat 
exchangers  would  yield  greater  performance  values  when  compared  to  their 
macro-scale  counterparts.  This  meant  that  much  higher  gradients  of  the 
concerning  parameters  were  expected  to  exist  within  the  solutions.  Therefore, 
the  need  for  a  suitable  grid  arose,  and  this  need  determined  how  fine  the  grid 
should  be.  It  was  a  significant  effort  to  obtain  valid  (grid  independent)  results.  The 
following  section  indicates  the  procedures  taken  to  achieve  grid-independent 
results. 

Two  prismatic  elements  were  considered  for  the  internal  part  of  the 
exchangers:  a  quadrilateral  element  and  a  triangular  element.  Following  the 
element  selection,  a  specific  measure  was  identified  to  properly  quantify  the 
quality  of  such  grids.  The  most  suitable  measure  for  this  task  was  the  pin 
perimeter  discretization.  Furthermore,  to  make  a  consistent  grid  comparison,  the 
grid  growing  factors  toward  the  planar  and  orthogonal  directions  to  the  pins’ 
cross  sectional  area  were  kept  constant.  Analogous  refinement  was  also 
conducted  within  the  entry  and  exit  sections  of  the  heat  exchangers. 

The  grid  independent  study  was  performed  with  the  X/D=1 .5,  S/D=1.5  and 
H/D=1  heat  exchanger  configuration  with  a  pin  diameter  of  100pm.  This 
configuration  was  selected  due  to  its  high  area  amplification  factor,  which 
indicates  a  very  high  pin  density.  It  was  anticipated  that  the  solutions  would  have 
high  gradients  in  all  concerning  parameters,  which  gives  the  solution  a  very 
strong  dependence  on  grid  refinement.  Thus,  to  achieve  grid  independence  the 
computed  heat  transfer  flux  had  to  reach  less  than  two  percent  variation  between 
subsequent  grid  refinements. 
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Table  4. 


Grid  independence  for  model  with  triangular  elements 


Parameter 

Grid  1 

Grid  2 

Grid  3 

Perimeter  division 

50 

70 

90 

#  of  cells 

88143 

141076 

201165 

#  of  nodes 

137380 

159516 

234304 

AQ  (W) 

0.00271 

0.0027 

0.00268 

Error  (%)  for  AQ 

1.12 

0.37 

0.00 

T  bulk  out  (K) 

305.49 

305.46 

305.44 

Table  5.  Grid  independence  for  model  with  quadrilateral  elements 


Parameter 
Perimeter  division 
#  of  cells 
#  of  nodes 
AQ  (W) 

Error  (%)  for  AQ 
T  bulk  out  (K) 


Grid  1 

Grid  2 

50 

70 

87450 

139802 

72030 

118800 

0.00257 

0.0026 

4.10 

2.24 

305.21 

305.32 

Grid  3 

Grid  4 

90 

110 

198660 

278018 

172284 

245504 

0.00265 

0.00268 

1.12 

0.00 

305.37 

305.43 

Tables  4  and  5  show  the  pin  perimeter  discretization  characteristics  and 
the  significant  data  extracted  from  the  simulations.  The  triangular  elements  were 
selected  because  they  proved  to  be  more  robust  and  because  the  solutions  with 
the  triangular  elements  required  less  grid  refinement  to  obtain  stable  outlet  bulk 
temperatures. 
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Figure  1 1 .  Comparison  of  heat  transfer  vs  grid  quality  for  the  examined 

models 

Figure  11  compares  the  results  from  simulations  performed  with  both,  the 
triangular  and  the  quadrilateral  element  models  versus  their  grid  quality  (number 
of  discretization  points  within  one  pin  perimenter.)  These  tests  were  performed 
for  the  case  of  Reynolds  number  of  100.  As  shown,  both  models  converged  to 
the  same  solution  of  0.00268  W.  However,  the  triangular  modeling  show  results 
that  overestimates  the  heat  transfer  when  using  a  grid  too  coarse  to  properly 
capture  the  high  gradients.  In  contrast,  the  quadrilateral  elements  underestimate 
the  heat  transfer.  Therefore,  selecting  the  triangular  elements  was  simple, 
specially  noting  that  grid  independent  results  were  obtained  with  20%  less 
discretization  points  around  the  pin  perimenter.  Furthermore,  the  deviation  from 
the  grid  independent  result  was  300%  larger  with  the  quadrilateral  elements. 

Figure  12,  plots  the  models’  performance  versus  the  time  required  for  the 
simulation  execution.  Selecting  the  number  of  divisions  per  pin  perimeter  was 
also  important.  From  both  plots  (figures  11  and  12,)  the  choice  of  70  points 
provided  the  best  trade-off  between  computational  time  and  accuracy  of  the 
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results.  At  70  points  per  pin  diameter,  the  results  predicted  a  heat  transfer  rate  of 
only  0.4%  different  than  the  finest  configuration  (using  110  points). 


Figure  12.  Comparison  of  the  outlet  temperature  achieved  vs.  computational 

time 


B.  THEORETICAL  VALIDATIONS 

Numerous  inspections  were  performed  in  order  to  validate  the  results 

generated  by  CFD  -  ACE.  The  contour  plots  of  temperature,  velocity  and 

pressure  were  examined  carefully  to  verify  that  they  satisfy  the  boundary 

conditions.  The  velocity  profiles  were  inspected  for  a  fully  developed  condition 

before  the  entrance  into  the  array.  In  addition,  the  output  file  generated  upon 

completion  of  every  simulation  was  carefully  examined  and  analyzed.  For 

example,  the  mass  flow  rate  computations  were  examined  for  an  exact 

agreement  with  the  inlet  and  outlet  boundaries.  This  was  verified  by  assuring  that 

the  imbalance  was  smaller  than  1%  of  the  total  value.  The  heat  transfer  summary 

was  also  carefully  inspected.  The  energy  flow  into  and  out  of  the  model  at  the 

corresponding  boundaries  was  examined  to  be  equal  to  the  total  amount  added 
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from  the  pins  and  the  endwall.  Also,  the  heat  imbalance  was  examined  to  check 
if  it  was  an  order  of  magnitude  smaller  than  the  net  energy  addition.  However, 
because  the  magnitude  of  the  non-dimensional  parameters  is  extremely  small 
(due  to  the  heat  exchangers’  dimensions,)  computation  errors  were  introduced. 

1.  Energy  Balance 

This  was  examined  for  agreement  between  the  bulk  outlet  temperature 
during  the  post-processing  of  the  results  and  the  temperature  that  was  computed 
from  a  simple  energy  balance  according  to  the  following  relation: 


Tout  Tn 


+  -^-  (22) 

mC„ 


The  AQ  was  computed  from  the  heat  transfer  summary  from  the  output  file 
runoff  each  simulation.  If  the  two  temperatures  had  more  than  a  0.1%  difference, 
the  simulation  was  repeated  with  a  finer  grid  or  with  an  increased  number  of 
iterations.  In  most  cases,  we  did  not  have  any  difficulty,  except  for  a  few 
simulations  at  a  Reynolds  number  of  1 ,000  and  for  some  of  the  configurations 
with  the  highest  area  magnification  factors. 

2.  Effectiveness 

The  independence  of  heat  exchanger  effectiveness  was  tested  by  varying 
the  prescribed  boundary  temperature  conditions.  Thus,  during  the  grid 
independence  tests,  three  different  boundary  temperatures,  306  K,  310  K,  and 
315  K  were  used  with  all  the  grid  refinements.  The  results,  shown  in  table  6,  were 
extremely  close  with  values  within  0.1%  from  each  other. 


Table  6.  Effectiveness  independence  of  boundaries’  temperatures 


Model 

1st  (50  points) 
2nd  (70  points) 
3rd  (90  points) 
4th  (110  points) 


Twall.l- 306 

0.868 

0.887 

0.896 

0.905 


Twall,2— 31 
0 

0.868 

0.887 

0.897 

0.905 


Twall,3—  315 

0.869 

0.888 

0.897 

0.905 


31 


For  grater  assurance,  heat  exchanger  No.  7  was  used  using  two  very 
distinct  boundary  temperatures,  306  K  and  400  K.  The  effectiveness  results  were 
again  within  0.1%  of  each  other. 

3.  Effectiveness  Based  on  NTU 

Effectiveness  is  associated  to  the  number  of  transfer  units  (as  shown  in 
equation  18).  Since  this  relation  is  valid  for  every  heat  exchanger,  it  can  be  used 
to  further  validate  the  computational  results.  In  order  to  thoroughly  test  the 
computational  model  against  the  trend  of  this  relation,  heat  exchanger  No.  7  was 
used.  Two  different  models  were  created,  one  with  two  rows  of  pins  and  the 
second  with  four  rows  of  pins..  The  computed  NTU  results  are  plotted  in  the 
purple  squares  of  figure  13. 


Figure  13.  Effectiveness  based  on  NTU  for  HX#7 
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The  black  line  in  Fig.  13  plots  the  computed  NTU  values  using  Eq.  18,  the 
results  show  how  close  the  CFD  models  approach  the  theory. 

4.  Compressibility  Effects 

Our  initial  ideas  of  simulating  the  problem  by  using  a  constant  density  for 
the  models,  similar  to  previous  numerical  works,  collapsed  very  early  since  some 
of  the  pin-fin  arrays  are  highly  packed.  It  was  anticipated  that  these  arrays  would 
require  high  differential  pressures  across  the  array,  which  breaks  down  the 
assumption  of  constant  density. 

Heat  exchanger  configuration  No.  7  was  used  to  simulate  various  mass 
flow  rates  with  and  without  variable  density  configurations.  The  tests  were 
performed  at  Re  =  100  and  1 ,000.  The  variable  density  model  used  the  ideal  gas 
law  to  compute  the  appropriate  density  values.  Figure  14  shows  these  results, 
which  show  significant  difference  in  their  friction  factors,  up  to  50%  difference  at 
Re  =  1,000. 


Consequentially,  variable  density  was  incorporated  in  all  subsequent  tests. 
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C.  LAMINAR  TO  TURBULENT  TRANSITION 

An  examination  of  the  transition  region  for  the  micro  heat  exchangers  was 
also  necessary.  For  this  case,  heat  exchanger  configuration  No. 9  was  used, 
because  the  results  can  be  corroborated  with  existing  experimental  data. 

1.  Turbulence  Modeling 

Within  our  computational  tool,  CFD-ACE,  various  turbulence  models  can 
be  used.  The  more  often  used,  the  “kt-£t”  model,  does  not  offer  very  accurate 
results  during  transition  because  it  only  considers  the  overall  effect  of  turbulence 
in  the  mean  flow  field.  However,  when  the  interest  occurs  with  parameters  such 
as  the  heat  transfer  coefficient  and  the  friction  factor  during  laminar-to-turbulent 
transition,  it  is  suggested  to  use  a  Low  Reynolds  number  model.  These  models 
permit  the  integration  of  momentum  and  “kt-£t”  equations  all  the  way  to  the  wall. 
The  difficulty  with  them  is  that  the  first  grid  point  must  be  placed  in  the  laminar 
sublayer  (y+~1 )  in  order  to  properly  connect  the  laminar  sublayer  with  the  log-law 
(wake)  layer.  Therefore,  this  kind  of  model  requires  the  use  of  very  fine  grids 
near  the  wall  boundaries. 

We  selected  the  Low  Reynolds  number  “kt-£t”  model  from  Chien.  For  more 
information  on  this  are  see  (Ref.  29)  or  (Ref.  30).  The  “kt-£t”  equations  have  been 
modified  to  include  the  molecular  viscosity  that  dominates  in  the  near  wall 
regions.  The  general  form  of  the  Low  Reynolds  number  models  is  given  by  the 
following  equations: 

|(pk)  +  /-(pu,k)  =  ^-((m  +  — )-f^-)  +  p(P-£-D)  (23) 

at  dx,  dx .  o,  5X; 

J  J  K  J 


at 


(p£)  +  —  (pUj£)  =  —  ((M  +  -^-)— J  +  f2^  +  E  (24) 

ax;  J  ax;  o'dx'  £i  1  k  E2  2  k 


Mt  =  C„f„  (25) 


The  Chien  model  parameters  appearing  in  the  equations  above  are: 


34 


(26) 


=  0.09  ,  CEi  =  1 .35 ,  C£2  =  1 .8 ,  CTk  =  1 .0 ,  ct£  =  1 .3 
fM  =1-exp(-0.0115y+),f1  =1.0 
Rp 

f2=1-0.22exp(-(— i)2) 

6 

OwU  C 

D  =  ^,E  =  -2v(4)exp(-0.5y+) 

y  y 

The  turbulent  kinetic  energy  was  expressed  using  the  turbulent  intensity 
variation,  which  is  provided  by  the  following  correlation  (Ref.  29): 

I  =  0.16Re'1/8  (27) 

for  the  core  of  a  fully  developed  duct  flow. 

As  suggested  in  (Ref.  29),  within  internal  flow  problems,  the  duct  hydraulic 
diameter  is  a  good  choice  for  a  turbulent  length  scale.  Further 
assumptions/settings  included  a)  zero  roughness  conditions  were  set,  b)  to  avoid 
numerical  instabilities,  the  upwind  first  order  scheme  was  used  for  advection,  and 
c)  the  number  of  iterations  of  the  solution  process  was  increased  in  order  to 
achieve  the  same  convergence  criteria  as  in  the  laminar  cases. 

2.  Transition  Region 

As  the  Reynolds  numbers  increased,  it  was  anticipated  for  the  turbulence 
parameters  to  increase.  Thus  in  each  simulation,  the  variation  of  the  y+  at  the 
wall  boundary  was  examined.  In  order  for  y+  to  remain  in  the  order  of  unity  further 
grid  refinement  was  required.  A  y+  <  1  was  achieved  at  the  lower  endwall 
boundary,  but  at  the  pins’  surface,  y+  increased  with  height  and  with  downstream 
location  from  the  stagnation  region  as  shown  in  Fig.  15.  The  larger  values  of  y+  at 
increasing  height  were  due  to  the  extrusion  power  law  that  was  used  for  the  Z 
dimension  during  the  grid  construction.  The  increase  in  y+  along  the  pins’ 
perimeter  can  be  explained  by  the  acceleration  of  the  flow  due  to  the  diffuser’s 
shape  region  between  pins,  thus  causing  a  local  increase  of  the  turbulence 
levels.  Figure  15  results  were  obtained  at  a  Reynolds  number  of  3,500. 
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Figure  15.  y+  variation  with  the  half  pin  perimeter  at  several  heights 


Despite  further  grid  refinements  and  increased  number  of  iterations,  it  was 
not  possible  to  reduce  the  y+  values.  Hence,  turbulent  models  require  very  fine 
grids,  a  lot  of  computer  memory,  and  even  greater  amount  of  computing  time. 
Figure  16  represents  a  three  dimensional  carpet  plot  of  y+  for  a  z-axis  cut  of  the 
model  for  a  Reynolds  number  of  3,500  at  z=0. 


Figure  16.  y+  variation  in  the  lowest  cells  layer  of  the  model 
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Also  very  encouraging  from  a  heat  transfer  point  of  view  is  that  the 
maximum  heat  fluxes  at  the  pins’  surface  are  met  at  the  35-40  degrees 
downstream  of  the  stagnation  point.  The  y+  in  this  region  is  still  increasing, 
because  it  is  the  initial  region  of  the  flow  acceleration.  As  expected  from  theory, 
very  similar  to  the  y+  trends  were  met  for  the  turbulent  kinetic  energy. 

We  try  to  identify  the  transition  in  a  plot  that  compares  the  effective  array 
heat  transfer  coefficient  with  the  required  specific  fluid  friction  power.  This  plot  fits 
our  purpose  well,  since  it  combines  both  heat  transfer  and  fluid  flow 
characteristics.  Figure  17  provides  this  information  for  both  the  laminar  and  the 
turbulent  models  for  the  Reynolds  numbers  between  1 ,000  and  8,200. 


h  n  vs  E  for  transition  region  in  HX-9 
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Figure  17.  Comparison  of  laminar  and  turbulent  models  for  HX-9 


Both  models  overlap  until  Re  >  2,000.  After  that  point  the  laminar  model 
under-predicts  heff.  In  the  configuration  of  Hamilton  (Ref.  7),  the  fall  of  the  laminar 

model  occurred  at  Re  >  1,000.  Therefore,  first  estimates  indicate  that  transition 
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occurs  in  different  regions  for  different  configurations  and  depends  on  the  area 
density  of  every  array.  For  heat  exchangers  with  smaller  area  magnification 
factors  the  transition  region  is  shifted  toward  higher  Reynolds  numbers.  At  even 
higher  Reynolds  numbers  shock  waves  are  formed.  This  will  be  examined  more 
extensively  in  the  following  chapter. 

The  Nusselt  number  follows  similar  trends  as  the  heat  transfer 
coefficients,  therefore  Nusselt  numbers  decay  near  and  after  the  shock  regions. 
The  shock  wave  problem  was  somewhat  expected  due  to  the  micro-scale 
character  of  the  exchangers.  When  comparing  micro-  with  macro-scale  heat 
exchangers  at  similar  Re  numbers,  the  micro-scale  heat  exchangers  have  to 
operate  at  higher  velocities,  which  give  rise  to  the  shock  problems.  But  since 
CFD-ACE  is  not  a  robust  compressible  solver,  it  is  not  very  reliable  for  analyzing 
problems  at  these  high  velocities.  Consequently,  we  use  the  results  of  the 
turbulent  model  until  a  Re  <  6,200,  where  the  maximum  Mach  number  was 
approximately  0.83. 

Another  aspect  that  has  to  be  examined  is  the  variation  of  the  effective 
viscosity.  We  found  that  the  effective  viscosity  reached  and  exceeded  values  five 
times  greater  than  the  laminar  viscosity  in  the  core  of  the  flow  (upper  symmetry 
plane  boundary  condition)  in  a  regime  of  the  Reynolds  numbers  around  2,000. 
Figure  18  presents  a  carpet  plot  with  the  effective  viscosity  values  for  a  z-axis  cut 
of  the  array  just  below  the  upper  symmetry  boundary  of  the  heat  exchanger  for 
the  Reynolds  number  around  2,000. 
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Figure  18.  Variation  of  effective  viscosity  in  the  HX-9  for  Re=2000 

It  is  generally  suggested  that  when  this  criterion  has  been  reached,  the 
modeling  must  switch  into  the  turbulent  regime. 

D.  EXPERIMENTAL  CORROBORATION 

1.  Validation  with  Historical  Data 

As  discussed  above,  the  current  CFD  results  were  compared  to  existing 
experimental  data.  Such  experiments  used  the  same  configuration 
(X/D=S/D=2.5,  H/D=1)  from,  Chyu  (Ref.  10)  and  Metzger  et  al.  (Ref.  6).  Chyu  et 
al  (Ref.  12)  conducted  an  experimental  work  based  on  a  well  known  mass 
transfer  technique,  the  naphthalene  sublimation  experiments,  and  through  the 
use  of  the  mass  transfer  analogy,  obtained  heat  transfer  results.  The  difference 
between  Chyu  and  the  current  study  are  the  following  definitions:  Chyu’s 
Reynolds  number  was  defined  as  follows: 


where  Vmax  is  the  average  flow  velocity  in  the  minimum  flow  area  (between 
a  pins  row)  and  D  is  the  actual  pin  diameter. 

Chyu  defined  the  Nusselt  number  as  follows: 
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where,  again,  D  is  the  pin  diameter  and  h  is  the  heat  transfer  coefficient 
that  comes  from: 


h  =  q"(Twall-Tbulk)(30) 

where  q"is  the  heat  flux  and  Tbuik  is  the  bulk  temperature  in  the  channel 
based  on  the  mass  transfer  analogy.  Because  it  was  impossible  to  equate  their 
Tbuik  using  the  information  provided  in  their  paper,  and  because  the  bulk  outlet 
temperature  used  in  the  current  study  (the  net  outlet  temperature  of  the  coolant 
fluid,)  was  not  consistent  with  theirs,  it  is  better  to  compare  their  Nu  number 
results  with  a  slightly  modified  version  of  our  Nu  number  (where  the  pin  diameter 
is  used  instead  of  the  hydraulic  diameter.)  The  results  can  be  compared  by 
examining  Table  7  and  Figure  19. 


Table  7.  Data  comparison  of  present  study  with  Chyu  et  al  (1999). 


Reph 

Nuoh 

1033 

19.7 

1498 

24.7 

2067 

30.0 

2790 

36.4 

3617 

42.7 

4650 

48.7 

6200 

51.8 

6613 

51.3 

8267 

36.9 

Present  study  Re 
equivalent 
969 
1406 
1939 
2617 
3393 
4362 
5816 
6204 
7754 


Present  study  Nu 
equivalent 

14.3 

17.9 
21.8 

26.4 

30.9 
35.3 

37.5 
37.2 
26.7 


Nu  from  Chyu  et 
al  correlation 

17.6 

21.9 

26.4 

31.5 

36.6 
42.4 
50.1 
52.0 
59.3 


The  correlation  that  Chyu  provides  is  Nu  =  0.32Re0  583  for  Re>6000.  The 
above  extracted  data  is  from  an  extension  of  their  correlation  toward  lower  Re 
numbers  because  our  micro  heat  exchanger  could  not  reach  their  higher 
Reynolds  due  to  the  shock  wave  issues. 
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Nu  vs  Re  comparison  with  data  from  Chyu  et  al  (1999) 


■  Present  study  Nu  equivalent  | 
Extented  Nu  correlation  from  Chyu  et  al 
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Figure  19.  Nu  vs.  Re  comparison  with  historical  data  for  HX-9 


If  we  consider  that  Chyu  tested  an  array  with  7  rows  of  pins  and  the 
current  CFD  model  for  HX-9  had  only  6,  the  difference  in  the  definitions  of  the  Nu 
number  and  other  parameters  (such  as  the  viscous  dissipation  not  being  included 
in  this  study,)  the  agreement  and  the  scale  of  the  findings  is  sufficiently 
satisfactory.  Moreover,  the  parallel  trend  that  the  two  curves  of  the  compared 
data  are  following  is  a  very  positive  result. 

By  comparing  our  findings  with  Metzger  et  al.  we  again  achieved  a  very 
similar  trend,  with  smaller  differences  of  about  15%.  Metzger  used  the  same 
definition  as  Chyu  et  al.  for  the  Re  and  Nu  numbers,  while  his  heat  transfer 
coefficient  was  derived  from  the  equation 


h  = 


q/A 


wetted 


(^wall  "  ^ref  ) 


(31) 


where  for  the  tref  definition  and  computation  the  reader  is  referred  to  (Ref. 
6).  The  fact  that  the  computational  model’s  trend-line  develops  in  the  region 
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between  the  two  predictions  of  the  two  empirical  researchers,  is  very  positive  for 
the  performance  fidelity  of  the  model. 


4fi 


Nu  vs  Re  comparison  with  Metzger  [1982)  data 
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Figure  20.  Comparison  of  Nu  as  a  function  of  Re  with  Metzger’s  results 


Also  Metzger  in  his  study,  found  that  the  average  Nusselt  number 
decrease  after  the  sixth  row  is  negligible  for  this  configuration  and  it  can  be  taken 
as  if  it  has  reached  a  steady  value  of  an  infinite  length  array. 

In  an  earlier  study  Chyu  (Ref  28),  also  researched  the  pressure  losses 
and  the  friction  factor  for  the  same  specific  configuration.  He  defined  the  friction 
factor  as  follows: 


2AP 
pV  2N 

r'  w  max 


(32) 


where  AP  is  the  pressure  difference  between  the  inlet  and  exit  sections  of 
the  array,  N  is  the  number  of  pin  rows  and  the  Vmax  is  defined  as  above.  His 
findings  are  presented  in  figure  20. 
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Figure  21.  Friction  factor  vs  Re  from  Chyu  [1989]  for  X/D=S/D=2.5,  H/D=1 


The  present  study  produced  equivalent  parameters  depicted  in  figure  21  - 
agreement  is  excellent!  Chyu’s  experimental  array  consisted  of  7  rows  of  pins 
and  appears  that  this  difference  did  not  have  the  influence  that  the  Nu  number 
results  showed. 


Present  study  f  variation  based  on  Chyus  definitions 


Friction  factor  variation  of  the  turbulent  model  for  HX-9 
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Figure  22. 


Marques  et  al.  (Ref.  16)  tested  the  same  configuration  as  Chyu  et  al.,  but 
in  the  micro-scale.  Marques’  pin  diameter  was  500pm,  the  same  as  to  the  current 
study.  Accordingly,  the  HX-9  design  is  identical  with  the  design  of  Marques,  but 
with  one  row  less.  Also  because  Marques  et  al.  heated  only  the  one  endwall,  we 
could  not  compare  the  present  study’s  heat  transfer  findings  with  theirs.  However 
a  comparison  of  the  friction  factor  could  be  very  useful.  After  they  corrected  the 
parameters,  they  validated  their  results  against  Chyu’s  study,  as  shown  in  figure 
23. 


Figure  23.  Friction  factor  as  a  function  of  Re  number  from  Marques  [2004] 


Consequently,  there  is  a  very  similar  trend  for  the  values  of  the  friction 
factor  predicted  above  which  seems  to  tend  towards  a  constant  value  at  Re  > 
6,000.  This  corresponds  to  a  value  of  Re  >  5,000  in  the  current  study. 

2.  Validation  with  NPS  Experimental  Research 
Roussakies  (Ref.  17)  and  Bosserman  (Ref.  18)  conducted  experimental 
research  using  identical  configurations  to  the  current  study.  Roussakies  worked 
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with  configurations  exactly  similar  to  HX-1,  2,  10  while  the  Bosserman  tested  all 
of  the  current  configurations  except  for  HX-7,8  and  9  due  to  problems  with 
blockage  effects  of  the  arrays’  passage  ways.  Therefore,  only  HX-1  and  10  were 
validated  against  their  results. 

Bosserman’s  results  of  the  friction  factor  for  some  of  the  HXs  like  3,  4  and 
6,  in  which  he  found  numbers  of  15  <  f  <  100  for  Re=100,  indicate  trends  that  can 
only  be  explained  by  the  heat  exchangers  being  clogged  (not  properly  released.) 

Figure  24  plots  the  Nu  for  the  first  heat  exchanger  for  a  range  of  Re. 


Comparison  of  Nu  vs  Re  for  HX-1 


Re 


Figure  24.  Comparison  of  Nu  vs.  Re  for  HX-1  with  experimental  work 


From  the  plotted  data,  we  have  good  agreement  for  Re  >  600,  where  the 
points  are  shown  converging  in  a  similar  trend.  We  noticed  a  similar  deviation 
with  the  Nu  results  for  the  average  array  heat  transfer  coefficient  when  it  was 
plotted  as  a  function  of  specific  fluid  friction  power.  This  was  expected,  since  both 
parameters  are  in  great  dependence. 
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Figure  25.  Comparison  of  harray,ave  vs  E  for  HX-1  with  experimental  work 

Regarding  the  friction  factor,  our  findings  were  again  corroborated.  Our 
results  were  found  to  be  between  the  trends  of  the  experimental  results.  Even 
more  positive  is  the  fact  that  the  data  showed  similar  decay  as  the  Re  increased 
-  shown  in  figure  26. 


Comparison  of  f  vs  Re  for  HX-1 


26 

25 


1  8 


s  15 

s 

13 

10 

oa 

0  5 
03 
00 


|  »  Present  compuational  model  | 

|  ~  Roussakies's  results  (2004)  | 

|  A  BotsarniaiVs  results  (2005)  ] 


— - - 1 - - - 1 - - - 1 - - - 1 - - - 1 - » 

100  200  300  400  500  flOO  700  800  800  1000  1100 


Re 


Figure  26.  Comparison  of  f  vs  Re  for  HX-1  with  experimental  work 
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The  comparison  of  HX-10  with  experimental  data  led  to  stronger 
corroboration,  especially  with  Roussakies’s  findings.  His  Nu  variation  with  Re  is 
closer  to  the  computational  model  results  than  Bosserman’s. 
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Comparison  of  Nu  vs  Re  for  HX-10 
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Figure  27.  Comparison  of  Nu  vs  Re  for  HX-10  with  empirical  data 


Roussakies  predicted  values  for  Re  >  400  close  to  the  numerical  model’s 
trend-line.  In  his  data,  the  same  behavior  is  seen  in  the  average  heat  transfer 
coefficient  values.  This  happens  according  to  figure  28,  for  low  friction  power 
values,  as  compared  to  Bosserman’s  results,  which  also  predict  a  well 
established  agreement  for  the  friction  factor.  The  trends  between  Roussakies 
and  our  model  in  figure  29  are  extremely  similar  and  the  values  are  very  close. 
The  predicted  friction  factor  values  are  almost  identical  for  a  region  of  Re  >  400. 
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Figure  28. 


Figure  29. 
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Consequently,  the  computational  models  of  the  heat  exchangers  that  have 
been  developed  and  are  discussed  in  this  chapter,  produced  results  very  similar 
to  earlier  experimental  data  in  the  macro-scale,  where  a  suitable  correction  was 
applied,  and  extremely  similar  where  the  accurate  equivalency  from  their 
definitions  was  known  and  has  been  implemented.  As  to  the  comparison  with  the 
micro-scale  results,  we  are  sufficiently  satisfied  with  the  existent  corroboration. 

Therefore,  based  on  the  validation  by  both  the  micro  and  the  macro  scale 
data,  the  results  of  the  macro-scale  experiments  can  be  directly  applied  to  micro¬ 
scale  heat  exchangers. 
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IV.  RESULTS-DISCUSSION 


Here  we  discuss  selected  performance  characteristics  of  the  heat 
exchangers  that  were  analyzed  (shown  in  Table  3)  which  for  convenience  are 
divided  into  two  sets  according  to  their  pin  diameter.  The  characteristics  are 
chosen  to  highlight  various  features  of  the  behavior  of  the  heat  exchangers. 

A.  HEAT  TRANSFER  FLUX,  RATE,  AND  COEFFICIENT 

The  heat  transfer  coefficient  is  one  of  the  most  critical  and  interesting  of 
the  examined  parameters.  Since  it  was  not  possible  to  directly  examine  the 
behavior  of  the  local  heat  transfer  coefficient  in  CFD  -  VIEW,  we  inspected  heat 
fluxes  and  heat  transfer  rates  instead. 


Pins  and  enwalls  performance  comparison  for  HX-3 


Figure  30.  Pins  and  endwall  heat  transfer  comparison  for  HX-3 


HX-3  was  found  to  be  one  of  the  best  performers,  since  it  had  a  very 
packed  configuration  but  a  bigger  pin  diameter.  The  contribution  of  the  endwall 
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and  the  pins  to  the  heat  transfer,  along  the  total  length  of  the  array  is  compared 
in  figure  30  for  Re=500.  It  can  be  seen  that  for  the  first  row  of  pins  the  heat 
transfer  rate  achieved  values  more  than  250  percent  greater  than  the  endwall, 
which  was  reduced  greatly  for  the  following  rows  without  any  great  alteration  in 
the  percentage  difference  contribution.  To  confirm  that  result,  we  should  inspect 
the  heat  flux  distribution  at  those  wall  boundaries.  The  endwall  showed  very  high 
fluxes  for  the  very  initial  part  of  the  Hx,  but  reduced  very  quickly  and  with  a  very 
high  gradient  along  the  stream-wise  distance  as  is  presented  at  the  figure  31 . 
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Figure  31 .  Endwall’s  heat  flux  distribution  in  HX-3 


After  the  first  row  the  values  of  the  actual  wall  flux  were  falling  to  a  level  of 
10,000  W/m2.  On  the  other  hand  when  we  examined  the  heat  flux  according  to 
pin  chord-length  in  the  first  row  for  several  heights,  we  noticed  that  it  reached 
near  maximum  values  from  almost  a  quarter  of  the  half  height  which  was 
modeled.  This  means  that  the  pin  keeps  the  highest  achieved  heat  flux  values  for 
more  than  75  percent  of  its  length.  Furthermore,  the  maximum  values  for  the 
heat  flux  reached  are  around  35-40  degrees  downstream  from  the  flow  direction 
angle  and  are  maintained  until  100  to  105  degrees,  from  which  point  they  start  to 
reduce  dramatically  due  to  generation  of  vortex  shedding  and  recirculation 
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effects.  Figure  32  is  a  plot  of  this  behavior  with  angle  from  the  forward  flow 
direction  to  360  degrees  around  a  pin. 


Heal  flux  distribution  across  pin  perimeter  for  Isr  row's  pin  in  HX-3 


Figure  32.  First  row’s  pin  heat  flux  distribution  for  several  heights  in  HX-3 


It  is  worth  noting  the  scale  of  the  heat  fluxes  presented  and  the  great 
potential  of  the  micro  heat  exchangers  to  dissipate  large  heat  fluxes.  The  values 
of  heat  flux  that  the  pins  reach  are  extremely  high  and  were  highest  for  the 
second  row  as  illustrated  in  figure  33. 
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Heat  flux  z-axis  profile  for  HX-3  In  35  downstream 


Figure  33.  Heat  flux  z-axis  profile  for  35°  downstream  the  pin  surface 


After  the  second  row,  the  values  of  the  heat  flux  were  reduced,  but  not 
with  rates  as  high  as  those  for  the  endwall.  The  second  row  was  the  one  that 
clearly  achieved  the  greater  values  and  this  is  probably  happen  due  to  the 
acceleration  of  the  flow  after  passing  the  first  row,  due  to  the  nozzle  effect  which 
drives  it  more  toward  the  next  row’s  pin  center.  In  figure  34,  which  presents  the 
contour  distribution  of  the  heat  flux  in  a  first  and  a  second  row  pins,  the  red 
region  of  higher  flux  values  has  shifted  toward  the  stagnation  point  for  the  second 
row  pin  compared  to  the  one  of  the  first  row.  However,  the  first  row  finally 
produced  bigger  heat  transfer,  when  we  look  back  at  the  comparison  plot  for  the 
pins  and  endwall.  This  could  be  explained  by  the  reason  that  the  first  row  pin 
achieved  greater  area  average  heat  flux  values. 
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Figure  34.  Heat  flux  contour  plot  of  1st  and  2nd  row  pins  in  HX-3 


By  averaging  the  heat  flux  values  at  these  two  pin  surfaces  using  the 
calculator  tool  of  CFD  -  VIEW  we  confirmed  that  result,  but  as  we  checked  and 
for  other  arrays,  we  found  that  it  did  not  happen  in  all  the  configurations.  Also, 
because  the  great  difference  in  the  heat  transfer  participation,  that  pins  and 
endwall  resulted  it  was  important  to  inspect  how  other  configurations  behaved.  In 
a  comparison  for  a  Reynolds  number  equal  to  500  for  HX-1,  we  got  opposite 
results  to  those  in  HX-3.  In  HX-1,  the  endwall  provided  more  heat  transfer,  but 
not  with  as  much  difference  as  the  pins  did  in  HX-3. 
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Pins  and  endwalls  performance  comparison  for  HX_1 
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Figure  35.  Pins  and  endwall  heat  transfer  comparison  for  HX-1 


The  greater  differentiation  is  found  in  the  contribution  of  the  first  row,  and 
as  we  checked  the  following  rows,  we  noticed  that  this  difference  dies  quickly, 
mostly  due  to  the  more  general  wall  behavior.  The  pins  continue  to  contribute  a 
high  enough  amount,  as  compared  to  their  initial  amount  of  heat  transfer.  At  this 
point,  we  should  mention  that  HX-1  did  not  return  the  greater  heat  transfer 
coefficient,  as  compared  to  the  other  arrays.  To  get  a  total  picture  of  this  matter, 
we  examined  HX-2  which  has  the  increased  span-wise  ratio  S/D.  The 
participation  in  this  array  was  comparable  for  both  -  the  endwall  produced  greater 
heat  transfer  in  the  first  row  only  while  the  pins  prevailed  in  the  following  rows 
with  their  more  constant  performance. 
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Pins  and  endwalls  performance  comparison  for  HX-2 


Figure  36.  Pins  and  endwall  performance  comparison  in  HX-2 


The  pins  produced  a  20  -  30  percent  greater  heat  transfer  from  row  2  to 
row  10.  Consequently,  the  pins  and  endwall  contribution  is  not  constant;  in  fact,  it 
is  strongly  dependent  upon  the  configuration.  More  specifically  if  we  examine  the 
area  magnification  factor  of  the  inspected  arrays,  we  can  state  with  great 
assurance  that  the  pins’  participation  depends  on  the  area  density  of  every  array. 
Thus,  in  HX-3,  with  the  higher  area  density,  the  pins  contribute  the  greatest  in  a 
comparison  of  the  three  arrays. 

Furthermore,  in  both  HX-1  and  HX-2,  the  second  row  of  pins  produced  the 
greater  heat  transfer  rate  in  contrast  to  HX-3.  Accordingly,  which  row  produces 
the  maximum  average  flux  and  subsequently  which  row  dissipates  more  heat  is 
something  varied,  and  depends  upon  the  configuration  too.  However,  the  heat 
flux  in  all  of  them  reached  its  maximum  local  values  in  the  second  row.  The  result 
that  the  second  row  achieved  the  higher  local  heat  fluxes  is  supported  by  figure 
37  which  is  a  plot  of  the  flux  distribution  at  the  upper  symmetry  plane  along  the 
pin  chord-length  for  all  rows  of  HX-1 . 
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Comparison  of  row’s  pin  heal  flux  distribution  for  HX-1 
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Figure  37.  Heat  flux  distribution  for  the  pins  of  HX-1  at  the  core  of  the  flow 


However,  in  order  to  make  a  more  fair  judgment  of  the  pins  and  endwall 
contribution,  we  use  a  chart  that  illustrates  the  attained  average  heat  fluxes  of 
every  row,  which  are  independent  of  the  flow  wetted  area  factor,  as  in  figure  38. 


Average  heat  flux  of  pins  and  endwall  for  HXsl.2,3 
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Figure  38.  Heat  flux  comparison  of  pins  and  endwall  in  HXsl ,  2,  3 
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Here  we  can  observe  that  the  difference  in  the  first  three  configurations 
varies  from  50  to  250  percent,  based  on  the  endwall  values.  Moreover,  it  shows 
more  clearly  that  the  performance  of  the  pins  is  more  generally  dependent  on  the 
heat  exchanger’s  wetted  area  magnification  factor.  In  addition,  it  is  clear  that  the 
second  row  does  not  always  accomplish  the  highest  flux  as  happens  in  HX-3. 

For  the  average  heat  transfer  coefficient  for  the  two  sets  of  arrays,  figures 
39  and  40  present  their  variation  with  Reynolds  number.  It  is  clear  that  HX-3  and 
7  achieve  the  greatest  values  in  the  higher  Re  number  regime  for  the  two  sets. 


h  vs  Re  for  HXs  1-4.9 


Figure  39.  Average  heat  transfer  coefficient  vs  Re  for  HXs  1-4,9 


From  a  more  total  view,  exchangers  6,  7,  and  8  reached  the  higher 
coefficients,  which  might  have  been  expected  since  they  have  the  greatest  area 
density.  When  we  inspected  the  first  set,  HX-4  (X/D=3,  S/D=1.25)  was  found  to 
have  greater  performance,  compared  to  HX-2  (X/D=1.25,  S/D=3),  suggesting 
that  the  stream-wise  spacing  variation  proved  to  have  a  greater  effect  compared 
to  span-wise  distance.  This  result  is  consistent  with  the  theory,  since  in  HX-4,  the 
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flow  is  forced  to  follow  a  more  circuitous  path  around  the  pins  and,  as  a  result,  to 
increase  the  interaction  with  them.  In  HX-2  with  its  wider  passages,  the  gain  in 
the  heat  transfer  is  not  so  significant,  since  the  flow  has  greater  direct  through- 
flow  regions  at  the  side  of  the  pins.  Also  worth  noting  are  the  values  for  HX-3 
which  are  far  higher  in  the  figure  than  all  the  others.  The  combined  reduction  in 
both  spacings,  stream-wise  and  span-wise,  was  the  main  thing  that  forced  the 
third  configuration  to  perform  so  much  better  than  the  others  and  to  attain  a 
similar  performance  as  HX-6  and  8  with  the  small  pin  diameter. 


vs  Ref<w  HXs  5-8,10 


Re 


Figure  40.  Average  heat  transfer  coefficient  vs  Re  for  HXs  5-8,10 


In  the  second  set  of  Hx,  the  pin  diameter  was  around  three  times  less, 
which  resulted  in  an  increased  H/D  ratio.  This  was  due  to  the  constant  height 
channel  that  we  kept  for  all  the  configurations.  In  comparing  the  two  plots  of  the 
two  sets,  it  is  easily  noticeable  that  the  main  trends  of  the  performance  that  have 
been  commented  on  previously  are  qualitatively  the  same  for  the  second  set,  but 
with  a  great  change  in  between  the  configurations’  quantitative  relation.  This  is 
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very  significant,  since  it  proves  that  the  alteration  of  one  of  the  three  ratios,  the 
pin  height  to  diameter  ratio,  had  a  significant  effect  on  the  resultant  performance. 
Thus  the  optimization  of  the  micro  heat  exchangers  is  not  so  simple  since  it 
involves  a  complicated  interaction  of  the  different  length  ratios  in  the  Hx. 

It  is  also  important  to  note  that  the  increase  in  the  H/D  ratio  had  an 
amplification  effect  on  all  of  them,  but  with  a  different  factor  for  every  one  of  the 
exchangers  as  seen  from  the  above  discussion. 

A  chart  plotting  the  heat  transfer  coefficient  versus  specific  fluid  friction 
power,  based  on  wetted  area  would  be  very  useful  in  order  to  evaluate  the 
performance. 


h  .  vs  E  for  HXs  1-4,9 


Figure  41 .  Performance  comparison  for  HXs  1  -4,9 


The  optimum  Hx  would  be  one  that  maximizes  the  heat  transfer  rate  while 
minimizing  the  frictional  losses  in  the  flow.  For  the  first  set  of  arrays  figure  41 
provides  that  information.  Instead  of  comparing  the  average  heat  transfer 
coefficient,  we  selected  the  effective  one,  because  it  better  represents  the  actual 
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heat  transfer  gain  from  every  configuration,  since  it  was  derived  by  taking  into 
account  the  magnification  due  to  the  array  area  density.  The  plotted  points  on 
each  line  are  for  increasing  Reynolds  number  values  of  100,  200,  300,  500,  700, 
and  1 ,000  in  the  direction  of  increasing  E. 

Between  HXs  1,  2  and  9  there  is  no  great  difference,  except  that  the 
second  one  which,  for  the  same  Re  values,  reaches  slightly  higher  values  of  heat 
transfer  coefficient,  with  the  analogous  required  energy,  than  the  others  would 
demand  operating  at  a  similar  point.  However,  HX-3  and  HX-4  have  a  great 
difference  between  them,  with  the  fourth  starting  from  a  value  for  the  heff  of  450 
W/m2/K  and  reaching  a  maximum  of  1,000  and  the  third  starting  from  that  point 
and  reaching  the  amount  of  4,000  W/m2/K  for  the  same  Reynolds  number.  If  we 
consider  that,  for  the  heat  transfer  coefficient  near  1,000,  one  demands  100 
times  more  power  to  drive  the  flow,  we  can  see  that  there  is  a  great  difference 
between  their  performance. 

HXs  2  and  4  are  closer  with  the  fourth  one  performing  better  until  a  value 
of  Re  equal  to  500,  whereas  after  this  point,  HX-2  produces  slightly  (maximum  30 
percent)  smaller  values  of  heat  transfer  coefficient  while  demanding  around  three 
to  four  times  less  energy  to  operate. 

In  the  second  set  of  inspected  configurations,  all  the  plotted  performances 
have  moved  toward  the  right  upper  region  of  the  chart  except  for  HX-10. 
Furthermore,  it  is  interesting  that  their  differentiations  are  scaled  similarly  in  both 
figure’s  directions  for  all  of  them  except  HX-7.  The  most  impressive  characteristic 
of  this  comparison  is  the  heat  transfer  coefficient  values  that  HX-7  reached 
considering  its  higher  area  density.  But  the  tremendous  amount  of  fluid  friction 
power  needs  to  be  compared  to  all  the  other  arrays.  The  deterioration  in  its 
performance  was  clearly  noticeable  when  it  approached  the  shock  wave 
operation  region,  as  the  heat  transfer  reduced  while  the  required  specific  power 
grew  from  200,000  to  almost  500,000  W/m2. 
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h„  vs  E  for  HXs  5-8.10 


Figure  42.  Performance  comparison  for  HXs  5-8,10 


To  increase  the  heat  transfer  performance,  it  is  much  more  cost  effective 
to  move  upward  in  the  diagram  to  a  more  packed  configuration,  instead  of 
changing  the  flow  attributes  and  increasing  the  Reynolds  number. 

Also,  if  we  consider  HX-1  and  10,  which  have  almost  identical 
performances  but  differ  a  lot  in  their  geometrical  characteristics,  the  gain  from 
increasing  the  H/D  ratio  in  the  tenth  has  been  diminished  by  the  increase  of  the 
other  two  spacings,  X/D  and  S/D,  as  compared  to  the  first  array.  This  is  an 
indication  that  the  height  to  diameter  relation  has  a  stronger  influence  compared 
to  the  other  two  ratios,  and  thus  needs  more  attention  in  the  optimization 
process. 

Another  way  to  inspect  the  heat  transfer  performance  of  the  tested 
configurations  is  to  compare  them  against  the  performance  of  the  similar 
dimensions  empty  duct.  Thus  in  this  case,  it  can  be  examined  the  difference  that 
the  presence  of  the  pin  made,  not  only  by  contributing  with  their  surface  but 
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enforcing  the  flow  to  follow  the  more  tortuous  paths  around  them  and  increase 
the  heat  dissipation  of  the  wall  boundaries. 
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Figure  43.  Variation  of  arrays’  effective  heat  transfer  coefficient  over  the  empty 

duct  for  HXs  1-4,9 


Figures  43  and  44  provide  the  ratio  of  the  effective  array  heat  transfer 
coefficient  over  the  corresponding  of  the  empty  duct  for  various  Reynolds 
numbers.  The  empty  duct  heat  transfer  coefficient  came  up  after  taking  in 
account  the  constant  Nusselt  number  of  the  infinite  width  square  duct,  which  is 
equal  to  7.54  for  the  thermally  developed  flow,  as  it  is  reported  by  Kays  and 
Crawford. 

The  configurations  like  HX-3,  6,  8  achieved  from  five  to  almost  thirty  times 
greater  heat  transfer  coefficients  when  consider  their  corresponding  empty  duct, 
but  the  HX-7  returned  values  of  the  order  of  thirty  to  more  or  less  one  hundred. 
This,  as  we  will  examine  the  friction  factor  in  a  following  section,  was  not  without 
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an  admirable  energy  requirement  that  almost  forbids  any  thoughts  of  applicability 
that  can  this  configuration  would  provide. 

Furthermore,  the  fact  that  for  the  arrays  the  heat  transfer  coefficient 
increases  significantly  by  increasing  the  flow  characteristics  (Reynolds  number), 
which  does  not  occur  so  intensively  for  the  case  of  the  empty  duct,  is  also 
important. 


Figure  44.  Variation  of  arrays’  effective  heat  transfer  coefficient  over  the  empty 

duct  for  HXs  5-8,10 


B.  NUSSELT  NUMBER  VARIATION 

The  variation  of  Nusselt  number  with  the  Reynolds  number  is  shown  in 
figure  45  for  the  first  set  of  arrays.  The  trends  are  not  very  different  compared  to 
the  behavior  of  the  average  heat  transfer  coefficient. 
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Nu  vs  Re  foi  HXs  1-4,9 


Figure  45.  Nusselt  number  comparison  for  HXs  1-4,9 


The  actual  relation  between  the  heat  exchangers  has  changed  a  little 
because  the  Nusselt  number  has  been  influenced  by  the  hydraulic  diameter,  as 
defined  in  a  previous  chapter,  according  to  equation  2.  This  differentiation  first 
concern  HX-2,  which  has  moved  in  a  lower  scale,  and  HX-9,  which  has  moved 
upward  relative  to  HX-1,  always  compared  with  the  heat  transfer  coefficient 
trends.  This  is  valid  with  the  variation  of  their  hydraulic  diameter,  which  for  HX-2 
with  the  smaller  volume  and  the  same  cross  sectional  area  as  HX-1,  results  in  a 
smaller  diameter.  For  HX-9  the  reduction  of  both  ratios  (X/D  and  S/D)  produced  a 
smaller  fractional  reduction  in  the  volume  than  in  the  cross  sectional  area,  and 
consequently,  the  enlargement  of  the  hydraulic  diameter.  Furthermore,  when 
considering  the  Nusselt  number  variation,  the  opening  between  the  HX-  3  and 
HX-4  has  been  decreased,  since  the  hydraulic  diameter  of  567pm  of  HX-4, 
compared  to  the  307pm  of  HX-3,  appears  as  an  amplification  factor  for  HX-3. 
This  trend  is  clearer  in  the  second  set,  where  the  increase  of  the  H/D  ratio 
increased  the  heat  transfer  coefficient  for  both  HXs  7  and  8,  but  the  further 

increase  to  the  hydraulic  diameter  for  HX-8  (Dh,8=354pm  vs  Dh, 7=1 37pm) 
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dominated  and  caused  not  only  the  gap’s  total  diminishing  but  also  the  eighth 
array’s  outperformance  over  the  seventh  until  an  Re  equal  to  300. 

By  continuing  with  the  exchangers  in  the  second  set,  as  was  done  for  the 
heat  transfer  coefficient,  the  performance  curves  in  figure  46  include  the  effect  of 
the  increase  in  the  H/D  relation  and  the  alteration  to  their  trends  by  our 
characteristic  length,  the  hydraulic  diameter. 


Nu  vs  Re  for  HXs  5-8.10 


Figure  46.  Nusselt  number  comparison  for  HXs  5-8, 1 0 


In  judging  all  the  configurations’  performance  considering  the  Nusselt 
number,  against  the  heat  transfer  performance  of  an  open  channel  with  a 
constant  surface  temperature  boundary  which  produces  an  average  Nu  equal  to 
7.54  (Ref  26),  almost  all  of  them  attained  much  higher  values  for  the  whole 
tested  range  of  the  Re  number.  The  only  exceptions  were  HX-1  and  HX-2,  which 
cross  the  limit  of  the  constant  open  channel  performance  around  a  Reynolds 
number  of  200. 
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C.  FRICTION  FACTOR  EFFECTS 

The  most  obvious  effect  that  occurred  in  the  friction  factor,  in  comparing 
the  two  sets,  was  the  magnification  due  to  the  increase  of  the  H/D  fraction  that 
resulted  from  the  reduction  of  the  pin  diameter,  while  keeping  the  channel  height 
constant.  The  reduction  of  the  pin  diameter  caused  a  decrease  in  the  hydraulic 
diameter  and  also  a  great  increase  in  the  average  array  velocity.  However,  these 
two  changes  were  not  enough  to  offset  the  decrease  of  total  array  length  and, 
more  significantly  the  increase  in  the  pressure  gradient  between  the  entrance 
and  the  exit  of  the  exchanger. 


Friction  factor  f  vs  Re  for  HXs  1-4.9 
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Figure  47.  Effect  of  Reynolds  number  in  friction  factor  of  HXs  1-4,9 


Furthermore,  for  the  friction  factor,  like  the  other  parameters,  we  note  that 
the  amplification  due  to  the  alteration  of  the  H/D  was  not  proportional  for  all  the 
arrays.  Thus,  this  is  another  proof  that  the  changes  to  the  geometric  fractions  do 
not  act  independently  and  must  be  considered  together. 
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If  we  consider  the  difference  in  the  effects  between  HXs-2  and  4,  or  6  and 
8,  since  they  have  the  same  hydraulic  diameter,  and  thus  the  same  average 
velocity,  the  main  cause  of  it,  is  the  great  difference  between  the  two  differential 
pressures  of  the  compared  arrays.  Moreover,  this  cause  is  consistent  with  the 
theory  given  that  the  smaller  S/D  ratio  for  HXs-4  and  8  produces  a  narrower 
span-wise  spacing  for  the  pins  and  consequently  the  flow  streamlines  negotiate  a 
more  serpentine  path  around  the  pins,  which  leads  to  the  greater  pressure  drops. 
On  the  other  hand  in  HXs-2  and  6  there’s  a  more  direct  path  and  more  smooth 
streamlines,  and  result  in  smaller  energy  amounts  being  needed  to  drive  the  flow 
through  their  interior. 


Friction  factor  f  vs  Re  for  HXs  5-8,10 


Figure  48.  Effect  of  Reynolds  number  in  friction  factor  of  HXs  5-8,10 


Another  aspect  worth  noting  is  the  high  friction  factor  that  HXs  4  and  8 
produce,  especially  as  compared  with  HXs-3  and  7.  The  main  reason  for  that  is 
the  dependence  on  the  squared  of  the  average  velocity  at  every  Re  number  and, 
the  increased  hydraulic  diameter. 
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By  comparing  the  pressure  drop  across  the  inlet  and  the  exit  of  every 
configuration,  with  the  pressure  drop  required  from  the  corresponding 
dimensions  empty  of  pins  duct,  to  run  the  flow  through  it,  we  get  figures  49  and 
50.  The  amount  of  energy  required  for  the  several  arrays  in  order  to  increase  the 
heat  transfer  performance  shows  similar  behavior  to  the  corresponding  friction 
factor. 


Figure  49.  Ratio  of  HXs  1-4,9  pressure  drop  over  the  empty  duct  case  vs.  Re 


When  considering  an  optimization  and  a  selection  process  of  the  most 
suitable  HX  for  an  application,  the  amount  of  the  pressure  drop  required  for 
returning  the  maximum  heat  transfer  performance  or  the  inverse,  has  to  be 
examined.  Thus,  HX-7  requires  around  five  times  more  pressure  difference  of  the 
increase  in  heat  transfer  that  provides,  HX-8  proved  to  require  six  times  more, 
while  HX-6  proved  to  be  more  economical  by  demanding  three  times  greater 
pressure  difference  than  the  times  that  increase  the  heat  transfer  when 
compared  to  the  open  duct  case. 
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Figure  50.  Ratio  of  HXs  5-8,10  pressure  drop  over  the  empty  duct  case  vs.  Re 


D.  MACH  NUMBER  VARIATION 

The  deterioration  of  the  heat  transfer  performance  of  HX-7  is  a  very 
interesting  manifestation  of  Mach  number  variation.  If  we  also  consider  the 
extremely  small  values  of  the  hydraulic  diameter  for  these  micro-scale  devices  it 
is  obvious  that  high  velocities  have  to  be  simulated  to  provide  the  desired  tested 
Reynolds  numbers. 

Recognizing  that  CFD  -  ACE  is  not  the  most  suitable  solver  for  studies 
that  involve  Mach  numbers  over  0.7-0. 8  the  presented  results  and  effects  are 
rough  estimates  and  must  be  seen  as  only  adequately  accurate  in  revealing  the 
main  trends.  A  compressible  solver  like  CFD  -  FASTRAN  would  be  required  to 
resolve  shock  behavior  with  greater  accuracy  but  was  not  used  here  since  it  was 
peripheral  to  the  main  study. 

However  by  using  the  CFD  -  VIEW,  we  tried  to  identify  the  maximum 
Mach  numbers  that  occurred  in  every  array  and  to  plot  them  with  respect  to  the 
Reynolds  number.  They  are  presented  in  figures  51  and  52.  Suitable  correlations 
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to  define  this  operational  region  would  be  useful  in  the  design  process.  The 
coefficients  of  a  power  law  correlation  below  are  shown  in  table  8. 


M  =  C1ReC2  (33) 


Maximum  Mach  vs  Re  for  HXs  1-4.9 


Figure  51 .  Maximum  Mach  number  variation  with  Re  number  for  HXs  1-4,9 


Maximum  Mach  vs  Re  for  HXs  5-8.10 


Figure  52.  Maximum  Mach  number  variation  with  Re  number  for  HXs  5-8,10 
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We  have  to  mention  again  that  these  correlations  are  not  intended  as  an 
accurate  prediction  of  the  flow  Mach  number,  but  rather  as  a  rough  guess  in 
order  to  assist  in  an  attempt  to  avoid  any  interference  with  shock  waves. 


Table  8.  Mach  number  correlations’  coefficients 


HX# 

Ci 

c2 

1 

0.0002 

0.9644 

2 

0.0002 

1.009 

3 

0.0006 

0.9445 

4 

0.0003 

1.0442 

5 

0.0002 

0.9635 

6 

0.0003 

1.0158 

7 

0.0005 

1.0179 

8 

0.0011 

0.991 

9 

0.0002 

0.968 

10 

0.0002 

0.9606 

HX-7  demonstrated  these  shock  conditions  at  the  high  inlet  velocities  that 
the  theory  dictated  for  the  tested  Re  numbers  and  found  to  have  maximum  Mach 
number  values  of  about  1.1.  Furthermore,  if  we  consider  the  inlet  velocity  of 
60m/sec,  or  the  average  of  121  m/sec,  or  the  maximum  value  of  the  fully 
developed  profile  of  90m/sec  for  the  corresponding  Re  equal  to  1 ,000,  and  the 
drastic  reduction  of  the  available  flow  area  between  the  pin  passageways,  shock 
waves  in  these  nozzles  like  sections  seem  unavoidable.  The  contours  of  the 
Mach  number  and  pressure  distribution  around  the  throat,  illustrated  in  figure  53, 
is  another  indication.  We  recognize  that  the  coarser  unstructured  grid  around  the 
centerline  is  not  appropriate  for  capturing  the  occurring  shock  conditions. 
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Figure  53.  Contour  plot  of  pressure  and  Mach  number  distribution  in  the  4th 

row  of  HX-7 


In  addition,  the  fact  that  the  critical  pressure  ratio  of  0.528  when  a  shock 
can  occur  is  reached  is  a  clear  indication.  For  all  these  reasons,  we  are  confident 
that  the  shock  behavior  was  captured  to  some  degree  by  ACE  and  is  reflected  in 
a  decline  in  the  heat  transfer  coefficient  and  in  Nusselt  number  values  around  a 
Reynolds  number  of  1 ,000. 

Bosserman  (Ref.  18)  observed  a  very  similar  condition  for  HX-3,  the 
denser  of  the  first  set  in  his  experimental  research.  Also  HX-3  contains  the  same 
area  percentage  reduction  between  the  pins  as  HX-7  does.  He  observed  similar 
decreased  performance  conditions  around  a  Reynolds  number  of  1 ,000  at  the 
same  parameters.  Figure  54  presents  an  indication  of  the  variation  of  Nusselt 
versus  Reynolds  number  for  HX-3  with  his  data. 
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Figure  54.  Nu  with  Re  variation  for  HX-3  [from  Bosserman,  2005] 


In  considering  the  data  in  table  9  that  provides  the  required  velocities,  inlet 
duct,  and  average  velocities,  we  concluded  that  a  shock  condition  must  have 
happened  at  a  smaller  Reynolds  number  than  expected,  since  at  a  Re  equal  to 
700,  the  theory  gives  an  inlet  velocity  of  17m/sec  for  which  it  is  very  difficult  for 
the  flow  to  be  chocked  under  the  geometrical  features  of  this  array.  However,  if 
we  consider  the  blockage  affects  that  Bosserman  faced  as  well,  it  is  very 
possible  that  the  flow  was  at  higher  velocities  and  that  he  experienced  shock 
conditions  that  are  consistent  with  the  gathered  data. 


Table  9.  Inlet  duct  and  average  array  velocities  for  HX-3  for  several  Re 


R©Dh,ar 

m dot  ar(kg/s) 

Uave ar(m/s) 

Vin(m/s) 

100 

1 .5463E-06 

5.39 

2.68 

200 

3.0926E-06 

10.78 

5.36 

300 

4.6389E-06 

16.18 

8.05 

500 

7.7315E-06 

26.96 

13.41 

700 

1 .0824E-05 

37.75 

18.77 

1000 

1 .5463E-05 

53.92 

26.82 
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E.  SLUG  VELOCITY  PROFILE  CONSIDERATIONS 

One  of  our  considerations  was  whether  the  performance  of  the  heat 
exchangers  would  change  and  if  so  how  much,  by  applying  a  more  realistic  and 
not  fully  developed  velocity  profile  as  an  input  boundary.  Therefore,  the  problem 
changed  to  a  thermal  constant  surface  entry  length  combined  with  a  flow 
developing  condition.  Also  it  is  well  known  (Ref.  26)  that  the  thermal  layer 
generally  develops  faster  than  the  velocity  boundary  layer  for  fluids  with  a  Pr  <  5. 

Thus,  we  chose  five  arrays  from  both  groups,  with  the  smaller  and  the 
bigger  pin  diameter,  in  order  to  check  the  effect  of  all  the  spacing  ratios,  testing 
them  in  the  new  conditions.  We  removed  the  entry  duct  from  the  model  and  we 
applied  directly  a  “slug”  constant  velocity  at  the  inlet  of  the  exchanger  that 
corresponded  to  the  desired  Reynolds  number.  It  is  very  clear  from  figure  55  that 
the  performance  did  not  change  too  much  when  the  slug  inlet  velocity  models 
were  compared  with  the  fully  developed  velocity  profile  models.  We  can  observe 
a  slight  increase  in  the  heat  transfer  coefficient  at  all  the  ranges  of  the  Reynolds 
number,  with  another  slight  increase  for  the  energy  required  to  drive  the  flow 
through  the  exchanger. 
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vs  E  for  comparison  of  HXs  2-4, 6,8 


Figure  55.  Fully  developed  and  slug  velocity  profile  comparison  effect  at  the 

HXs  2-4, 6, 8  for  heff  and  E 


The  resulting  behavior  is  consistent  with  the  fluid  physics,  since  from  the 
flow  perspective  at  the  initial  part  of  the  endwall,  the  friction  gradient  will  be  very 
high  until  the  flow  will  be  developed  relative  to  the  endwall,  resulting  in  increased 
fluid  friction  power  and  also  in  an  enhanced  heat  transfer  performance. 

As  far  as  the  Nusselt  number  is  concerned,  it  is  higher,  especially  for  the 
lower  x+  values.  Also,  as  the  x+  increases,  the  Nusselt  number  has  to  converge 
to  the  values  of  the  thermal  length  only  problem,  since  the  effects  of  the 
combination  of  the  two  developing  conditions  are  diminished.  Thermal 
development  is  attained  when  the  x+  coordinate  reaches  values  of  0.1 . 


Table  10.  x+  values  for  HX  2-4, 6, 8  in  Re=100,1000 


HX#  Dh  ar(m)  L(m) 

2  0.000567  0.00625 

3  0.000307  0.00250 

4  0.000567  0.00600 


x+  for  Re=  1 00  x+  for  Re=  1 000 
0.312  0.031 

0.230  0.023 

0.299  0.030 
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6 

8 


0.000354  0.00208 

0.000354  0.00200 


0.166 

0.160 


0.017 

0.016 


In  table  10  we  notice  that  for  all  the  tested  exchangers,  the  thermal 
development  was  reached  for  Re=100,  however  we  are  slightly  below  this  point 
for  the  higher  Reynolds  number  flow  conditions.  Thus,  we  expected  the  deviation 
in  the  Nusselt  number  between  the  similar  geometrically  models  compared,  to 
slightly  increase  as  the  Reynolds  number  was  increasing,  which  did  not  occur. 


Nu  vs  Re  for  Comparison  of  HXs  2-4 .6,8 


Figure  56.  Nusselt  number  variation  for  compared  velocity  profiles  in  HXs  2- 

4,6,8 


The  most  probable  explanation  is  that  the  performance  characteristics  had 
changed  so  little  that  it  was  of  the  order  of  the  simulation  accuracy,  and  we  could 
not  identify  the  above  mentioned  characteristic  behavior.  Also  we  noticed  that  the 
selection  of  the  total  number  of  rows  for  modeling  every  configuration  was 
consistent  with  the  corresponding  x+  values  in  order  to  reach  a  near  thermal 
developed  condition,  as  was  explained  in  the  chapter  two. 
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A  plot  of  friction  factor  with  Re  in  figure  57  helps  to  compare  the  role  of 
inlet  velocity  profile  on  pressure  drop.  It  is  clear  that  there  is  very  little  influence  of 
the  velocity  profile  and  noticeable  only  for  HX-2  and  HX-6. 


f  vs  Re  tor  comparison  ol  HXs  2-4.6, 8 


Figure  57.  Friction  factor  variation  for  compared  velocity  profiles  in  HXs  2- 

4,6,8 


The  most  logical  explanation  is  that  for  the  heat  exchangers  that  provide  a 
larger  cross  sectional  area  for  the  flow,  the  effects  of  the  change  in  the  inlet 
velocity  profile  were  more  important  due  to  greater  interactions  with  the  endwall. 
In  the  other  configurations,  probably  the  presence  of  the  pins  and  their  dominant 
“denser”  arrangement  did  not  allow  great  changes  in  the  compared  parameters. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  objectives  of  this  study  were  successfully  met  by  computationally 
simulating  the  performance  of  ten  micro  heat  exchangers  with  varied  X/D,  S/D, 
and  H/D  ratios.  The  code  CFD  -  ACE  was  found  to  be  very  robust  and  user- 
friendly  and  was  used  with  great  confidence  for  the  subsequent  computations  of 
the  laminar  flow  regime  and  the  performance  maps  of  the  interesting  parameters. 
Most  of  the  significant  conclusions  have  been  pointed  out  during  the  discussion 
in  each  section/chapter. 

In  summary,  we  can  say  that  the  laminar  modeling  was  shown  to  be 
accurate  until  a  Reynolds  number  of  1,000.  At  higher  values  of  Re  there’s  a 
likelihood  of  transition  to  turbulence,  although  the  transition  value  is  not  constant 
and  is  greatly  dependent  on  the  geometric  configuration  of  each  array. 

The  ninth  and  tenth  configuration  was  well  enough  validated  by  historical 
macro-scale  and  current  experimental  micro-scale  data.  Thus,  we  can  conclude 
that  direct  application  of  the  macro-scale  dimensionless  parameters  can  be 
performed  in  the  micro-scale  dimensional  heat  exchangers  in  order  to  examine 
the  performance  of  any  desired  configuration. 

Moreover,  we  noticed  that  the  effects  of  an  alteration  of  any  of  the 
geometric  ratios  influence  the  resulting  effects  of  the  variation  of  the  other  ratios. 
This  was  more  pronounced  for  the  H/D  ratio. 

The  performance  maps  of  the  interesting  parameters  like  the  Nusselt 
number,  average  and  effective  heat  transfer  coefficient,  specific  fluid  friction 
power,  friction  factor,  and  Mach  number  were  also  provided  for  in  all  the 
configurations.  They  can  be  very  useful  tools  during  a  design  process.  The 
tremendous  ability  and  potential  for  heat  transfer  removal,  especially  for  the  more 
packed  arrays,  has  been  pointed  out  and  also  the  extreme  fluid  power  that  some 
of  them  required. 
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An  examination  of  the  inner  part  of  the  compact  heat  exchanger  was 
conducted,  and,  mostly,  the  behavior  of  the  pins  compared  to  the  endwall  was 
performed,  which  proved  the  great  contribution  that  the  pins  provide  in  the  total 
heat  transfer  removal  process. 

The  desired  tested  Reynolds  number  simulations  in  the  micro-scale 
dimensions  required  much  higher  velocities  which  led  to  shock  wave  conditions 
very  early  for  HX-7,  Thus,  the  computation  with  varied  properties  was  proved,  the 
best  selection.  Also  the  presence  of  vacuum  effects  at  the  exit  of  some 
exchangers  further  enhanced  more  the  choice  of  the  variable  properties 
modeling. 

The  alteration  of  the  inlet  velocity  profile  that  was  performed  by  a  constant 
“slug”  velocity  profile,  in  order  to  inspect  the  change  in  the  parameters’  trends, 
did  not  return  any  considerable  deviation  of  the  initial  maps. 

For  all  these  reasons,  the  optimization  and  the  design  of  micro-scale 
compact  heat  exchangers  is  a  complicated  process.  It  must  be  executed  with 
considerable  attention,  since  one  has  to  consider  more  physical  phenomena  than 
in  macro-scale  design. 

Likewise  to  get  a  more  total  picture  of  the  performance  in  future  work  we 
recommended  that  the  viscous  dissipation  effects  be  taken  into  account  in  a 
computational  model,  in  order  to  check  how  this  negatively  affects  the  output. 

Also  the  examination  of  configurations  with  other  shape  fins  that  will 
provide  better  flow  characteristics  and  will  be  more  effective  is  emerging,  while 
the  difference  on  the  heat  transfer  performance  that  the  liquids  will  provide  has  to 
be  considered. 

Finally  of  a  lesser  concern  is  the  fact  that  in  some  of  the  exchangers 
tested  (only  two  to  three)  for  the  higher  Reynolds  number  (700  and  more  for 
1,000),  the  residuals  had  a  periodic  pattern.  This  feature  is  an  indication  of 
unsteadiness,  which  we  recommend  should  be  checked  by  running  an  unsteady 
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simulation  to  inspect  any  effects  on  the  final  performance  behavior  of  the  heat 
exchangers. 
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